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Chapter 1 

The Coronal Heating Problem 

The Sun, our nearest star, is both the source of energy for hfe on Earth and a 
unique physics laboratory. Energy produced in the core of the Sun is transported to its 
surface and atmosphere. Prom here, through radiation and a flow of particles (the solar 
wind), it affects the Earth and forms the so-called heliosphere, a cavity inside the local 
interstellar medium which extends beyond the solar system boundaries. 

The visible surface of the Sun, the photosphere, has a temperature of about 6000 K 
and emits mostly visible light. High energy radiation (EUV, FUV, up to X-rays) is produced 
in the upper layers of the atmosphere: the chromosphere, the transition region and the 
Corona. On the other hand the steady production of this kind of radiation requires high 
temperatures, and the highest temperatures are localized in the Corona. Hence this gives 
rise to the so-called Coronal heating problem: Corona is counterintuitively much hotter 
(millions of degrees K) than the underlying photosphere. While there is general agreement 
that the magnetic field plays a fundamental role, and that the source of this energy derives 
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from convective motions in the photosphere (which have more than enough energy), the 
debate currently focuses on what physical mechanisms can transfer, store, and dissipate 
this energy between the photosphere and the corona. 

1.1 The Sun: its Interior and Atmosphere 

The solar interior is conveniently separated into zones by the different physical 
processes at work. The pressure, density, and temperature decrease going outward, from 
the center up to the surface, located at 1 solar radius (1 Rq ~ 7 x 10^ km). Energy is 
generated by nuclear reactions in the core, the innermost 25% of the radius. This energy 
is transported outward by radiation through the radiative zone and by convective flows 
through the convection region, the outermost 30%. 

The convection zone is the outermost layer of the solar interior. It extends from a 
depth of about 2 x 10^ km up to the visible surface. At the base of the convection zone the 
temperature is about 2 x 10^ K. This temperature is low enough that the heavier ions are 
not totally ionized. This increases the cross section for the radiation, ultimately trapping 
more heat, which in turns makes the plasma unstable to convection. Convection occurs 
when the temperature gradient required to radiatively transport the energy is larger than 
the adiabatic gradient, i.e. the rate at which the temperature would fall if a volume of 
material expanded adiabatically while rising toward the surface. Where this occurs a parcel 
of plasma displaced upward will be warmer than its surroundings and will continue to rise, 
i.e. convection instability sets in. Convective motions are quite efficient in carrying heat to 
the solar surface, and while the plasma rises it expands and cools. At the visible surface 
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Figure 1.1: Large field-of-view image, in visible light, of sunspots in Active Region 10030 
observed on 15 July 2002 by the Swedish 1-m Solar Telescope on the island of La Palma, 
Spain. It shows part of a sunspot group near disk center. Distance between tick marks is 
1000 km. Credits: the Royal Swedish Academy of Sciences. 

the temperature drops to 5800 K and the density to 2 x 10~^ gcm~^. 

The photosphere is the visible surface of the Sun. Since the Sun is made of ionized 
gas, this surface is not sharply defined, and it is actually a layer about 100 km thick (very 
thin compared to the 7 x 10^ /cm radius of the Sun). A number of features, which are 
relevant for coronal heating models, are observed in the photosphere, including the dark 
sunspots, the bright faculae, and granules. 

Sunspots (Fig. II. ip appear as comparatively isolated dark regions on the surface. 
They typically last for several days, although very large ones may live for several weeks. 
Sunspots are magnetic regions on the Sun characterized by a magnetic field of the order 
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of about 100 — 1000 gauss, noticeably stronger that the surrounding field (a few gauss). 
Sunspots usually come in groups with pairs of spots, with opposite field polarity. The field 
is strongest in the umbra, the darkest part of the sunspots, while it is weaker and more 
horizontal in the penumbra, the brighter part surrounding the umbra. Temperatures in 
the dark centers of sunspots drop to about 3700 K, compared to 5800 K for the surround- 
ing photosphere, due to the presence of the strong magnetic field which partially inhibits 
convective motions. 

Faculae are bright areas that are most easily seen near the limb, or edge, of the 
solar disk. These are also magnetic regions but the field is concentrated in much smaller 
bundles than in sunspots. While the sunspots tend to make the Sun look darker, the 
faculae make it look brighter. During a sunspot cycle the faculae actually win out over 
the sunspots and make the Sun appear slightly (about 0.1%) brighter at sunspot maximum 
that at sunspot minimum. 

Granules (Fig. 1 1.2] and also Fig. II. ip are small (their typical size is about 1000 km) 
cellular features that cover the entire Sun except for those areas covered by sunspots. These 
are the tops of convection cells where hot plasma rises up from the interior in the bright 
areas, spreads out across the surface, cools and then sinks inward along the dark lanes. 
Individual granules last for only about 8 minutes. The granulation pattern is continually 
changing as old granules are pushed aside by newly emerging ones (see movie from the 
Swedish 1-meter Solar Telescope). The average speed of the fiow within the granules is 
1 kms~^, but it can reach supersonic speeds of more than 7 kms^^ . This fiow, and its 
interaction with the magnetic field, generate waves which then propagates toward the upper 
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Figure 1.2: Image of a solar active region taken on 13 May 2003 near the center of the 
solar disk at hehographic coordinates NIO E15 degrees. The tick marks are 1000 km apart. 
The image is a filtergram taken in 430 nm "G-band" hght at the Swedish 1-meter Solar 
Telescope. The image was taken by Dr. Tom Berger of the Lockheed Martin Solar and 
Astrophysics Lab, Palo Alto, California. 

layers of the atmosphere. 

Supergranules are another flow pattern present in the photosphere, at a larger 
scale (about 35.000 km across) than granules. These features also cover the entire Sun 
and are continually evolving. Individual supergranules last for a day or two and have flow 
speeds of about 0.5 kms~^. 

The chromosphere is an irregular layer, 2000 — 3000 km thick, above the photo- 
sphere where the temperature rises from 4400 K (the temperature minimum) at its base 
to about 2 X 10^ K at the top. At these higher temperatures hydrogen emits light that 
gives off a reddish color (Ha emission). This colorful (chromo = color) emission can be 
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seen in prominences that project above the Hmb of the sun during total solar eclipses. The 
chromosphere is the site of activity as well. Solar flares, prominence and filament eruptions, 
and the flow of material in post-flare loops can all be observed over the course of just a few 
minutes. 

The transition region is a thin and very irregular layer of the Sun's atmosphere 
that separates the chromosphere from the much hotter corona. The temperature changes 
rapidly from 2 x 10^ -fC up to 1 x 10^ K over a 30 km distance. Hydrogen is ionized at 
these temperatures and is therefore difficult to detect. Instead the light emitted by the 
transition region is dominated by such ions as C IV, O IV, and Si IV that emit light in the 
far ultraviolet region (< 2000 A) of the solar spectrum that is only accessible from space. 

1.2 The Solar Corona 

The Corona is the Sun's outer atmosphere. Most of the visible light coming from 
the corona is light coming from the photosphere diffused trough Thompson scattering. 
This results in a faint emission compared to the photosphere, and it only becomes visible 
during total eclipses of the Sun, as shown in Fig. 11.31 or using a special instrument called 
coronagraph. The corona displays a variety of features including streamers, plumes, and 
loops. These features change from eclipse to eclipse and the overall shape of the corona 
changes with the sunspot cycle. However, during the few fleeting minutes of totality few, if 
any, changes are seen in these coronal features. 

Early observations by Harkness and Young in 1869 of the visible spectrum of the 
corona revealed bright emission lines at wavelengths that did not correspond to any known 
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Figure 1.3: Image of the solar corona in visible light taken on 29 March 2006 during total 
solar eclipse in Side, Turkey. Light from the photosphere is diffused through Thompson 
scattering. The image was taken by Koenraad van Gorp. 

materials. This led astronomers to propose the existence of "coronium" as the principal gas 
in the corona. The true nature of the corona remained a mystery until 1939, when Grotrian 
and Edlen identified the Fe XIV and Ni XVI lines in the corona. The corona is in fact 
heated to temperatures greater than 1 x 10^ K. At these high temperatures both hydrogen 
and helium (the two dominant elements), and even minor elements like carbon, nitrogen, 
and oxygen are totally ionized. Only the heavier trace elements like iron and calcium are 
able to retain a few of their electrons. It is emission from these highly ionized elements that 
produces the spectral emission lines that were so mysterious to early astronomers. 

The corona shines brightly in x-rays because of its high temperature. On the other 
hand, the cooler solar photosphere emits very few x-rays. This allows us to view the corona 
across the disk of the Sun when we observe the Sun in X-rays. To do this we must first 
design optics that can image x-rays and then we must get above the Earth's atmosphere. 
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which shields the Earth from this high energy radiations. In the early 1970s Skylab carried 
an x-ray telescope that revealed coronal holes and coronal bright points for the first time 
(these features were actually visible in earlier sounding rocket data, but they were not 
recognized). In the 1990s Yohkoh provided a wealth of information and images on the solar 
corona. Today SOHO and TRACE satellites are obtaining new and exciting observations 
of the Sun's corona, its features, and its dynamic character. 

Coronal loops are found above sunspots and active regions. These structures are 
associated with the closed field lines that connect magnetic regions on the solar surface. 
Many coronal loops last for days or weeks but most change quite rapidly. 

Coronal holes, where the corona is dark, were discovered when X-ray telescopes 
were first flown above the earth's atmosphere to reveal the structure of the corona across 
the solar disc. They are associated with "open" magnetic field lines and are often found at 
the Sun's poles (a region that is currently studied by the Ulysses mission). The high-speed 
solar wind is known to originate in coronal holes. 

Different layers of the solar atmosphere are characterized by different features, but 
most of these features are connected one another by the magnetic field. At solar minimum 
the solar magnetic field is approximately a dipole field, with the magnetic axis aligned to 
the rotational axis, so that we have "open" magnetic field around polar regions and closed 
magnetic structures around the equator. This picture is, of course, an approximation and 
the major departures are at the large scales the presence of an heliospheric current sheet (a 
region where the magnetic field changes rapidly its polarity), and at smaller scales a more 
disordered structure, which becomes more important during solar maxima. 
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The origin of the solar magnetic field is an active field of research. Not being a 
solid body the surface of the Sun is characterized by differential rotation, i.e. the poles 
have a lower rotation rate than the equator. In fact, the rotation period is of 35 days at 
the poles and 25 days at the equator. Helioseismology has demonstrated that differential 
rotation extends down to the convection zone. The deeper zones, the radiative zone and 
the core, appear to rotate as a solid body. The mechanism, called magnetic dynamo, which 
leads to the formation of the solar magnetic field seems to be connected with convection. 
This intense magnetic field rises toward the solar surface due to magnetic buoyancy and 
emerges at the photosphere. 

The organization of the magnetic field at the photospheric level gives rise to two 
different kind of regions, so-called "active regions" and "quiet-Sun regions", as shown in 
Fig. 11.41 All the solar surface is characterized by a magnetic field with a mixed polarity, but 
in active regions the unipolar areas are bigger and the magnetic field stronger. This structure 
is commonly called a magnetic carpet. There is general agreement that the active region field 
is the direct result of the large-scale solar dynamo. The continuous presence of quiet-sun 
regions also during solar minima, when the number of active regions considerably decreases, 
suggests that this flux might be generated by local dynamo action just below the sun's 
surface, driven by granular and supergranular flows [551 IISl EH EI] The opposite polarity 
regions of the magnetic carpet are connected by closed magnetic field lines, extending up to 
the upper layers of the solar atmosphere, called loops. In correspondence of active regions, 
where the magnetic field is stronger, these loops shine bright in EUV and x-rays as shown 
in Fig. 11.51 
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Figure 1.4: Le/t: Full disk magnetogram showing the strenght and polarity of the line-of- 
sight component of the photospheric magnetic field. The red rectangle highlights an active 
region, while the blue one highlights a quiet-Sun region. Right: Close-ups of the active and 
quiet-Sun regions. These magnetograms have been obtained with the Michelson Doppler 
Imager (MDI) onboard SOHO. 

The intimate relationship of the magnetic field and the structures of the outer at- 
mosphere is illustrated in Figure [L6l that shows images of the Sun at different wavelengths, 
corresponding to different temperatures and layers. In particular the blue continuum im- 
age ll.6al shows the photosphere, image ll.6bl is a magnetogram showing the line-of-sight 
component of the magnetic field (white for out-going directed field, black for in-going field), 
the upper photosphere shown in the light of Ca II K in image II. 6c^ the chromosphere in 
Ha line in image ll.6dl the transition region in the light of He II at 304 A in image II. 6e^ 
and then the corona seen in lines Fe IX 171 A formed at ~ 7.5 x 10^ K in ll.Gfl Fe XII 195 
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Figure 1.5: Left: Active regions 9628 and 9632 are rotating over the southwest hmb of the 
Sun in this image, showing the loops arching above them. This image was taken in the Fe 
IX/X 171 A channel at 23:59 UT on January 10, 2001. North is to the left, west to the top. 
The field of view is 512 arc sec N/S and 384 arc sec E/W. Right: This image of coronal 
loops over the eastern limb of the Sun was taken in the 171 A pass band, characteristic of 
plasma at 1 MK, on November 6, 1999, at 02:30 UT. Both images have been obtained by 
the Transition Region And Coronal Explorer (TRACE) satellite. 



A at ~ 1.5 X 10^ K in [L6gl Fe XV 284 A at ~ 2.5 x 10^ K in fOhl and 3 - 5 x 10^ K 
in the soft x-ray region in ll.6i[ The main information we have from these pictures is the 
correspondence between the regions of strong magnetic field at the photospheric level ()1.6bp 
and the activity in all the upper layers of the solar atmosphere up to the X-ray corona. The 
main structures seen in Figure [T.6al are the sunspots. Emergence of the magnetic field at the 
photosphere II. 6bl occurs in and near the spots and elsewhere. The upper photosphere 1 1 . 6c\ 
chromosphere [L6dl and transition region [T.6el show local brightening, heating, at the loca- 
tions of strong magnetic field. The coronal images [l.6fm.6il show a complex of loops. Note 
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(a) Blue Continuum (b) Magnetogram (c) Ca II K 




(g) Fe XII 195 A (h) Fe XV 284 A (i) Soft X-rays 

Figure 1.6: Multi- wavelength images of the Sun taken on 8 February 2001, all within hours of 
each other. Credits: ESA/NASA SOHO (a, b, e, h), TRACE (f, g), the Japanese-American 
YOHKOH satellite (i), and the Big Bear Solar Observatory (c, d). 
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that all space is not covered by loops and that hotter temperature loops tend to be nested 
inside loops of cooler temperatures. 

1.3 Overview of Coronal Heating Models 

In the 1930s Edlen, Grotrian, and Lyot established that the solar corona has a 
temperature of the order of 10^ K. After this, the fundamental step was done in the 1940s by 
Biermann [4] , Alfven [l] , and Schwarzschild [83] , who pointed out that the high temperatures 
of the corona are a direct consequence of the convective motions at the photospheric level. 
The convection does a mechanical work, which subsequently is transported and dissipated 
in the chromosphere and the corona. 

As pointed out in the previous paragraphs, the solar corona consists of magneti- 
cally confined regions characterized by closed structures (loops) roughly located around the 
equator, and regions of "open" magnetic field (coronal holes) roughly located around the 
northern and southern poles. Observations of the solar corona in EUV and X-rays show 
that great part of this emission takes place in the closed structures. This is the reason 
for which the "open" magnetic regions, which look faint in this high-energy range of the 
spectrum, are called coronal holes. In fact they look like "holes" in an otherwise bright 
corona. 

Whether the mechanism which heats open and closed magnetic regions may be 
similar or not, in this work we focus our attention on the magnetically confined regions of 
the solar corona. 

The active X-ray corona has a temperature of 1 — 5 x 10^ K, an electron and proton 
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numerical density of about lO^'^ cm~^, and is magnetically confined in the ~ 10^ — 10^ gauss 
magnetic field of active regions. These high temperatures are maintained by a heat input 
of about 10^ erg cm~'^ [95| . 

In the 1940s Biermann, Alfven, and Schwarzschild supposed that waves, generated 
at the photospheric level by convective motions, and then propagating upward into the 
corona, would have dissipated their energy leading to the high temperatures observed. 
Although the mechanisms which are able to transfer, store and dissipate the energy are still 
a matter of debate, the basic idea remains unchanged. The waves which are generated at the 
photospheric level include sound waves, gravitational waves, and magnetohydrodynamics 
waves. More recently it has been shown [64^ [871 l90l F78] that all but Alfven waves are 
dissipated and/or refracted before reaching the corona. Then, while these other waves 
contribute to chromospheric heating, it is only Alfven waves which are able to reach the 
corona. 

The Alfven wave is a purely magnetohydrodynamic phenomenon, and it is es- 
sentially an oscillation due to magnetic field line tension. In fact transverse motions of the 
magnetic field lines cause a force that tries to restore them to straight-line form. Linearizing 
the equations of magnetohydrodynamics (hereafter MHD) in the simple case of a homoge- 
neous plasma embedded in a homogeneous magnetic field Bq, Alfven waves are found as a 
transverse incompressible wave, propagating in the direction of the wave vector k with the 
dispersion relation: 

CO = v_Ak cos 9 (1.1) 
where oj is the wave frequency, f_4 = Bq/^/Attp the Alfven velocity, k the wave vector 
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modulus and is the angle between the magnetic field Bq and the wave vector k. Alfven 
waves due to photospheric motions are expected to have periods comparable to the 300 
seconds characteristic time of granules, whose characteristic dimension / is of the order of 
1000 km and their velocity is of the order of 1 km s~^. Shorter periods have been detected, 
but at noticeably reduced power levels. 

Coronal loop length is of the order of 10^ — 10^ km, with a typical sound speed of 
the order of 2 x 10'' cm s~^, and Alfven speed roughly 10 times larger, 2 x 10® cm s~^. An 
important parameter in plasma physics is the ratio between kinetic and magnetic pressure 
/3 = 8ttp/B^, that in the case of a coronal loop is of the order of /3 ~ 2 x 10 , i.e. a coronal 
loop is a magnetic dominated system. 

The basic problem with wave heating is that the wavelengths of an Alfven wave 
with a period typical of photospheric motions are too large to match coronal loop lengths. 
In fact from the dispersion relation (jl.ip (with v_a = 2x 10^ km s^^), for a period T ~ 10^ s 
it follows a wavelength A ~ 2 x 10^ km, which is of the same order as the length of the 
longest loops. Then they are quasi-static displacements of the magnetic fields rather than 
waves. 

Traditionally the mechanisms responsible for coronal heating have been divided 
into two main groups: AC (alternate current) and DC (direct current). 

1.3.1 High Frequency Models 

AC heating models propose two different mechanisms, Phase mixing (Heyvaerts &; 
Priest [43j) and resonant absorption (Davila [15]), in order to facilitate dissipation of these 
waves within about one wave-period. They both occur when Alfven velocity is not uniform. 
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Phase mixing occurs when Alfven waves propagate at different phase velocities 
along nearby fieldlines, making them come out of phase. 

Resonant absorption occurs whenever the Alfven velocity is nonuniform (e.g. if 
density is not uniform) in a loop cross section. The wave amplitude is enhanced in a narrow 
layer where the local Alfven resonance frequency matches the frequency of the global loop 
oscillations. Gradients in the magnetic and velocity fields are very large in this layer, and 
the wave energy is easily dissipated by Ohmic and viscous processes. 

1.3.2 Low Frequency Models 

Alternatively it has been proposed (Parker [68], [69], [70], [71]) that the X-ray 
corona is heated by dissipation at the many small current sheets forming in a coronal loop 
as a consequence of the continuous shuffling and intermixing of the footpoints of the field 
in the photospheric convection. The formation of these current sheets is conjectured by 
Parker [72] in the following way. Consider a region < x,y < I, < z < L, embedded 
in a uniform magnetic field aligned along the z-direction Bq = 3^6^ ■ Top and bottom 
boundaries are located at the planes z = and z = L. Supposing that in the plane z = 
we have a zero velocity field, and that in z = L an incompressible 2-D, i.e. Vz = 0, velocity 
pattern. A continuous mapping of the footpoints velocity pattern in the perpendicular mag- 
netic field {bx, by) is produced. This field spontaneously produces tangential discontinuities: 
the discontinuities appear in the initially continuous field at the boundaries between local 
regions of different winding patterns. The tangential discontinuities (current sheets) become 
increasingly severe with the continuing winding and interweaving, eventually producing in- 
tense magnetic dissipation in association with magnetic reconnection. Parker suggested 
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that this dissipation is largely in the form of bursts of rapid reconnection. It is this spo- 
radic explosive dissipation at the tangential discontinuities in the bipolar fields on the sun 
that creates the active X-ray corona. The heating occurs in bursts, which are estimated 
to involve individually 10^^ — 10^^ ergs. Such a burst is too small to be observable and he 
refers to the individual burst as a "nanoflare", because it is 9 orders of magnitude smaller 
than a large flare of 10^^ ergs cm— 2 

He eventually computes the energy input. Magnetic field lines connect the fixed 
footpoints at z = with the moving ones in z = L, which move about with the velocity 
pattern imposed. The field lines have more or less a uniform deviation 0(i) to the vertical, 
where 

vf 

tane(t)~y^ (1.2) 
ij 

supposing Q{t) < 1. The vertical component of the field is Bq, indicating with b± the 
orthogonal component we have 

6i = 5otane(0~ ^ (1.3) 

the field line tension opposes this movement with a stress of the order of 6_Li?o/4vr, so that 
the work for unitary time and surface (the power) done by the photosphere is 

... ^^^±^0 B^vH _2 -1 n A\ 

Vv ~ = — ^ — ergs cm s . (1.4) 

The input fiux then increases linearly with time. We know from observations that the input 
fiux is of the order of ~ 10^ ergs cm^"^ s^^. For a magnetic field = 10^ G, a velocity 
field V = Ikms^^, and a loop length L = 10^ km, it follows from equation (jl.4p that the 
observed input flux is reached at t ~ 5 x 10^ s, when b± ~ Bq/A and ~ 14° (the so-called 
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"Parker angle"). At this point it is conjectured that bursty rapid reconnection dissipates 
b± as rapidly as it is produced by the velocity forcing at the photosphere. In this way the 
input energy flux is on the average always of the order of ~ 10^ ergs cm~^ s"^. 
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Chapter 2 

The Reduced MHD Model for 
Coronal Heating 

A coronal loop (Figure [2. ip is a closed magnetic structure threaded by a strong ax- 
ial field, with the footpoints rooted in the photosphere. This makes it a strongly anisotropic 
system; the measure of this anisotropy is given by the relative magnitude of the Alfvenic 
velocity ~ 1000 kms'^ compared to the typical photospheric velocity Uph ~ 1 kms~^. 
So the photospheric velocity, that is the amplitude of the Alfven waves that are launched 
into the corona, is very small compared to the axial Alfven wave velocity. 

To investigate the properties and dynamics of a so complex system such as a 
coronal loop, we make a simplified model which captures the essential features of a real 
loop. The main (indeed defining) feature of a loop is its strong axial field which serves as a 
guide for the Alfven waves excited by photospheric motions. It is just the dynamics of these 
waves, propagating along the guide field, that we want to study. We first assume a simplified 
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Figure 2.1: Left: Image of coronal loops over the eastern limb of the Sun in the 171 A pass 
band (~ 1 MK) taken on November 6, 1999, at 02:30 UT by the TRACE satellite. Right: 
We use a simplified straightened out cartesian geometry to model a coronal loop. Top and 
bottom surfaces represent the two photospheric sections at which the loop is anchored. The 
plasma is embedded in a uniform and homogeneous axial magnetic field Bq. 

geometry, neglecting any curvature effect, and model a coronal loop as a "straightened out" 
cartesian box (Figure [2.ip . i.e. as a parallelepiped with an orthogonal square cross section 
of size i±, and an axial length L embedded in an axial homogeneous uniform magnetic 
field Bq. For a quantitative numerical study, we next adopt the so-called "reduced MHD" 
equations to model the dynamics of the plasma (Kadomtsev &: Pogutse [26], Strauss |88) 
and Montgomery [58ll59j ). Magnetohydrodynamics (MHD) is used to study long-scale low- 
frequency phenomena in plasma physics. When a plasma is embedded in a strong magnetic 
field, a further simplified set of equations ("reduced") are derived from the full set of MHD 
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equations, to model the dynamics of the plasma. 

In § 12.11 reduced MHD equations will be described in detail. In § 12.21 we give a 
brief review of current anisotropic MHD turbulence theory, a phenomenon naturally arising 
in an anisotropic environment such as a coronal loop threaded by a strong axial magnetic 
field. At last, in §[3]we describe the numerical parallel code that we have developed to solve 
the reduced MHD equations and that we have used to perform our numerical simulations. 



2.1 Reduced Magnetohydrodynamics 



The equations of reduced MHD have been derived in two different research fields. 
Kadomtsev & Pogutse [i6j and Strauss [88] have derived these equations in the context of 
fusion research to model the dynamics of a plasma embedded in a strong magnetic field. 
They were specifically thinking about tokamaks, but their derivation is not strictly limited 
to the geometry of these devices. Montgomery [58\ [59] has derived, in a different way, the 
same set of equations to study MHD turbulence in the presence of a strong magnetic field. 

The equations derived by all these authors are exactly the same. This is a clear 
indication that when a plasma is embedded in a strong axial magnetic field the turbulence 
which naturally develops is strongly affected by the magnetic field, acquiring its anisotropic 
features (Montgomery [581 [59], Shebalin et al. [H]), and that the overall dynamics is well 
described by the equations of reduced MHD (Kadomtsev & Pogutse |46j, Strauss [88]). We 
now derive these equations following Montgomery [59| . 
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The equations of incompressible resistive MHD are: 



Po{-^ + u-Vuj=-Vp+^^-^ + uV^u, (2.1) 
dB 

— = Vx{uxB) + ri—V^B, (2.2) 

V-u = 0, (2.3) 

V • S = 0, (2.4) 

j = X (2.5) 



where po is a constant and uniform mass density, u is the velocity field, B the magnetic 
field, j the electric current, p the kinetic pressure, c the speed of light. Using a simplified 
diffusion model, both the magnetic resistivity (rj) and the shear viscosity (z/) are constant 
and uniform. 

We now derive the set of equations for a plasma embedded in a strong magnetic 
field directed along the axial direction z. Let's suppose the following ordering for the 
magnetic field: 

B ^^Bo + + e S(^) + 0(e2), e < 1, (2.6) 

where 

Bo = Boe,. (2.7) 



We also suppose that gradient scales along the axial field direction are longer than perpen- 
dicular scales ~ e, an hypothesis which is actually one of the main results of modern 
anisotropic turbulence theory (see § I2.2.4p . We introduce similar expansions for velocity. 
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electric current and kinetic pressure: 

M ~ + £^(1) +0(e2), (2.8) 
j ^ j(0) +0{e'), (2.9) 

p ~ + ep(i) +0(e2). (2.10) 

Using as variables t, x, y, and C = ^z, the gradient operator can be written as 

V = e^d^ + Bydy + eezdc^ = V_l + ee^9^. (2.11) 

The magnetic field B is expressed in terms of the vector potential A : 

B = -^Boe, + V X A, (2.12) 

with 

A~ + eAW + 0(e2), (2.13) 

and we choose the Coulomb gauge V • A = for its algebraic convenience. 

Introducing the expansion (12.12P in the equation ()2.2p for the magnetic field, after 
a few algebraic manipulations we obtain: 

dA fl \ 

— = ux i-Boe^ + V X Aj -ijcj + Vx, (2.14) 

where x is a scalar potential, resulting from pulling off a curl operator from eq. (|2.2|) . for 
which we assume an expansion of the form 

X~^X0 + X^°^ + eX^'^+O(62). (2.15) 



The momentum equation (|2.ip becomes: 



Po{ ^ + u-Vu] = -Vp + - 5o - X + ^ X (V X A) + uV^u. (2.16) 
at J e c c 
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Considering the leading order contributions 0(l/e) from (|2.14p . (j2.15p . (|2.16p . and 0(1) 
from ([23|), ([23]) we have 

X e,So = -Vxxo, (2.17) 

jWxe,Bo = 0, (2.18) 

V^-ti(°)=0, (2.19) 

_^-(0) = ^V^xS(0) = -— ViA^. (2.20) 

Setting xo = Bq if, (12.17P gives 

uf'^ = V±(p X e^. (2.21) 
Equation (j2.18p imphes that the j has only the component along the axial direction 

j(°^=j.e,. (2.22) 
This with equation (j2.20p implies that A'^''^ also has vanishing orthogonal components 

A(°)=A^e2. (2.23) 
Consequently from (|2.12p for B^^^ we have 

= {V±A,) X e„ (2.24) 

so that B^*') has a vanishing axial component Bi^^ = 0. 

Considering the 0{1) terms in (j2.14p we have 

^ e, = uf> xe,Bo + uf x [{V±A,) x e,] - rjcj, e,+e,Bo^ + V^x^°^ • (2.25) 

Introducing the vorticity = V x u, and following the above expansion conventions, we 
have 

cj(o) = X u|l°^ = -e^ V'iif = u e^. (2.26) 
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After some algebra the 0(1) terms of equation (j2.16p can be written as 



Po 



, Mo) yy 1 



ez = l^o^ez + ^ (ef • Vx) e, + uVloo e, (2.27) 



Finally, considering the z-components of equations (I2.25P and (12.270 . we obtain the reduced 
MHD equations: 



PO 



(2.28) 
(2.29) 



Finally, from 0{e) of equation ()2.19p we have 



• + ^ = 0, (2.30) 



and considering the orthogonal component of equation ()2.25p 



X e^Bo = -V^x^^\ i.e. u^^ = ( ^ J x e,. (2.31) 



From this last equation it follows • u||^^ = 0, which inserted in (I2.30p implies that, in 
the absence of a parallel flow at the boundaries, we have 



ni°) = 0. 



(2.32) 
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2.1.1 Dimensionless Form and Boundary Conditions 

The equations of reduced MHD (|2.28p - (|2.29p written for the transverse fields u± 
and b_L in dimensional form are: 



/^^"^^r ^ ^ T7 f , , (b±-V±)b^ , Bpdb^ , 2 ,„oo^ 

Po [-g^ + (u^ ■ V^) J = - Vx + ^ J + ^ + + -^^-^ (2.33) 

^ = (b± • Vx) - • V^) + i?o^ + V^'^lb± (2.34) 

V_L • it_L = (2.35) 

V_L • b± = (2.36) 



where po is a constant and uniform mass density. To render the previous equations in di- 
mensionless form, we notice that the magnetic fields b± and Bq can be expressed in velocity 



dimensions by dividing by y/AirpQ, i.e. considering the Alfvenic velocities. As character- 
istic quantities we use the perpendicular length of the computational box i±, the typical 
photospheric velocity Up^, and the related crossing time t± = i±/uph- The dimensionless 
equations are then given by: 

+ {u^ ■ Vx) = - (p+^]+ {b± ■ V^) b^ + VA^ + ^^lu^ (2.37) 



dt ' " " V 2 y ^ " ^ ^ dz n 

^ = {b±- Vx) - {u^ ■ Vx) b^ + v_^^ + :^Vlb^ (2.38) 
dt dz TZm 

• u_L = (2.39) 

• 6^ = (2.40) 

where Vj^ is the ratio between the axial Alfvenic velocity and the photospheric velocity, i.e. 
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and 

1 z/ 1 rj(? 

are the kinetic and magnetic Reynolds numbers. 

Introducing the velocity and magnetic potentials and ijj 

= V X e,) 6^ = V X {i, e,) (2.43) 

w = (V X u_l)3 = -Vi(/9 (2.44) 
J = (V X 6^)3 = -ViV (2.45) 
equations (I2.37p - (l2.38j) in terms of potentials are written as: 

^ = ^^^ + [v^,V'] + :^ViV', (2.46) 

Ot OZ TZm 

% = ^y^% + b'' ^\ - ^\ + l^i'^' (2-47) 

where the poisson bracket of two functions / and g is defined as: 

[/,^] = ^^_^^. (2.48) 
' dx dy dx dy 

Using equations (I2.43P it can be shown that: 

[V',^'] = -(u±- V^)^, [j,^'] = (b±- Vx)i, [a;,</.] = (u^- Vx)w. (2.49) 
Using these relations, equations (|2.46p - (j2.47p can be rewritten as: 

^ + {u±-'^±)i^ = VA^ + ViV', (2.50) 

ot OZ Tim 

^ + (u± • Vx) a; = (bx • Vx) j + VA^^+ ^^l^^ (2-51) 
which are the dimensionless version of equations (|2.28p - (|2.29p . 
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In turbulence theory the fundamental variables are the Elsasser variables 

z^ = u_i±b_i, (2.52) 

in terms of which, and supposing the Reynolds numbers to have the same values TZ^ = TZ, 
the reduced MHD equations (p^Tjl - ipIiO]) are: 

^ = - (.- . V^) z+ + + 1 Vi^+ - V^P (2.53) 

^ = - (z+ . V^) z- - vj^ + 1 Vi.- - V^P (2.54) 
V_L • = (2.55) 

where P = p + h\/2 is the total pressure, and is linked to the nonlinear terms by incom- 
pressibility (j2.55p : 

2 

ViP = - 5^ (d,zj) {d,zt) (2.56) 

An analysis of equations (j2.53p - (|2.55p gives us a qualitative preview of the results 
which will be obtained both numerically and analytically in the following chapters. The 
linear terms in equations (|2.53p - (|2.55p 

dz"^ dz^ 

show that z^ fields present an Alfven wave propagation along the axial field direction. In 
particular z~ describes waves propagating in the direction of Sq, and z+ in the opposite 
direction; both at the Alfven wave velocity v_a. This wave propagation is present also 
when the nonlinear terms become important, and transport energy from the photospheric 
boundaries into the system. 

As boundary conditions at the photospheric surfaces {z = 0, L) we impose a ve- 
locity pattern which mimics photospheric motions. In terms of the Elsasser variables the 
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velocity is the sum 

2+ + z"=2tt_L. (2.58) 

In terms of characteristics this gives rise to a "reflection" of the Alfven waves at the bound- 
aries, where we can only impose a condition on the incoming wave (alternately z'^ and z^\ 
To mimic photospheric motions we impose a velocity pattern on the top and bottom planes. 
In terms of the Elsasser variables to impose a velocity pattern means using the constraint 

z-^^z- = 2uph (2.59) 

Since z"*" and z~ are, respectively, waves propagating toward the inside and the outside of 
the computational box, this is a "reflection" condition on these waves, i.e. 

z" = -z+ + 2uo at 2 = (2.60) 

z+ = -z~ + 2uL atz = L (2.61) 

At the boundary the value of the incoming wave is equal to the negative value of the 
outgoing wave plus twice the value of the velocity at the photosphere. 

A fundamental feature of the nonlinear terms [z^ ■ V^) z^ (and also the pressure 
term (j2.56p ) is the absence of self-coupling, i.e. the nonlinear interaction depends by the 
counter-propagating fields 2;^, and if one of the two fields were zero there would be no 
nonlinear dynamics at all. This is the basis of the so-called Alfven effect, which is described 
in § 12. 2. 2^ and is the basis of anisotropic turbulence phenomenology. 

2.1.2 Conservation Laws 

MHD theory provides a number of conservation laws which play an important role 
in turbulence theory. We are mainly interested in the Energy, Cross Helicity, and Magnetic 
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Helicity. Turbulence is usually studied with the hypothesis of periodicity in all three spatial 
directions. In our case along the direction of the axial field this condition breaks down. 
The flux terms which are usually neglected become important. In this paragraph we write 
the conservation laws for the three aforementioned quantities, including the flux terms and 
restricting our attention to the reduced MHD equations. 



Multiplying the momentum equation (j2.37p by u_l and the magnetic fleld equa- 
tion ()2.38p by we obtain: 



d_ 

dt 



hi 



(u^ ■V)u^-V±(p + 



{h± • V)u^ - (u^ • V)bj 



VA- 



V)b, + 



dz n 
1 



OZ TZ.m. 



(2.62) 
(2.63) 



On the other hand from 



B = + bj 



u = u_L, V • B = V • b_L = V • = 0, 



(2.64) 



it follows that 



(V X B) X B = 
(V X u_l) X 



-VB^ + (B • V) B = --Vb^ 

2 ^ ^ 2 ^ 



1. 



- Vui + (U_L • V) U_L 



du I 



V X (u^ X B) = (B • V) - (u^ • V) B = va^ + 

OZ 



(2.65) 
(2.66) 



+ (b^ • V) - (u^ • V) b^ (2.67) 
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In this way equations (|2.62p - (j2.63p can be rewritten as 



d (\ 



5^ 2"^'="^- 



5t U ^ 



b 



(V X ux) X - ( p + -ui 



\ 

V X (ux X B) + — Vibx 



+ (V X B) X B 
1 



7^ 



(2.68) 
(2.69) 



The first term between square brackets in the right hand side of equation (j2.68p is orthogonal 
to u_|_ and hence its scalar product with u_|_ is zero. Furthermore 



V • [B X (ui_ X B)] = (u_L X B) • (V X B) - B • [V X (u^ x B)] (2.70) 



and 



(u_L X B) • (V X B) = -u_L • [(V X B) X B] 



(2.71) 



B • [V X (u_L X B)] = -b^ • [V X (u^ X B)] 



(2.72) 



We can now write for the energy density, summing equations (I2.68P and (I2.69P 



— -ui + -bi 

at V 2 -^ 2 ^ 



-V • [B X (u_L X B) 



1 9 1 9 

u± • Vlu_L + — b_L • Vibj 



7^ 



7^ 



(2.73) 



The following relations hold: 



u± • Vj 



U_L • V^Ui^ ~ U_L • V^U_L = V • (u_L X Cj) — CJ^ 



bx • vib^ ~ bx • v^bi = V • (bx X j) - r 



(2.74) 

(2.75) 
(2.76) 



where 



cj = V X J = V X bj 



(2.77) 
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We can now write 



B X (u_L X B) + ( p + ^u5_^ - ^ (u_L X CJ + b_L X j) 

= -i{-^'+j') (2-78) 



The Poynting Flux S is given by: 



S = B X (u X B) = B^u - (B u) B 



(2.79) 



So using equation (|2.64|) we have 



S = {v\ + hi) ux - {h± ■ ux) {VA + b^) 



(2.80) 



and the component along the axial direction is 



5^ = S • = -VA (b_L • u_l) 



(2.81) 



At last we can write 



|(iui + ibi)+v.s.-i(.nf) 



(2.82) 



Calling the total Energy E, the Ohmic dissipation rate J, and the viscous dissipation rate 
J7 (which is the enstrophy divided by the Reynolds number TZ) 



E = \ J^d'x (ui + bi) , 



(2.83) 
(2.84) 
(2.85) 



and S the integral of the Poynting flux given by 



S* = - / d^x V • S = <i)daS • n = / da (u_l • b_L) - / da (u_l • b^) , (2.86) 

'v J Jz=L Jz=0 
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where the signs have been chosen so that S is positive when we have energy entering the 
system and negative when it leaves. We can write the integral of equation (I2.82p 

dE 



dt 



s-{n + j) 



(2.87) 



Another important conserved quantity is cross helicity H'^ , defined as: 



H' 



c 



d^x u_|_ • bj 



(2.88) 



V 



Multiplying the momentum equation (|2.37p by b_|_ and the magnetic field equation (j2.38p 
by u_L, after a few algebraic manipulations we obtain 



dt 



(u± • bj 



X (ux X bx) + (p+ ^ui jb^ + ^( ui + bi ]vAe, 



1 



+ ;^ ( b_L • Viu_L + u_L • Vib_L 



(2.89) 



When taking the integral over the volume the first and second terms in square brackets do 
not contribute because their only components lie in the orthogonal plane where we have 
periodic boundary conditions. But the third term cannot be neglected and we thus obtain 

^^J^d^xu^-h^ = ^ Jda {ul + hl){x,y,z = L)-{ul + hl){x,y,z = 0) 

+ 1 l^d^x (h^ ■ Viux + • Vibx) (2.90) 

Given a magnetic field B and its vector potential A for which B = V x A, in 
MHD the magnetic helicity is usually defined as 



H 



M 



d^x A • B, 



(2.91) 



but this definition is, in general, not gauge invariant. In fact taking a gauge transformation 
A' = A + Vx gives 



d X B • Vx = <J) da X B • n 

V Js 



(2.92) 
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Now equation (|2.9ip is gauge invariant, but only when the surface integral (j2.92p vanishes. 
This condition is satisfied when the normal component of the magnetic field vanishes at the 
boundary surface. For our coronal loop this condition does not apply, and definition (|2.9f |) 
cannot be used. An alternative expression has been proposed by Finn and Antonsen [25] 
(see also Berger and Field [3]) 

H^t, = [ d'x (A + Ao) • (B - Bo) (2.93) 
Jv 

where Bq = V x Aq is a reference field to be chosen suitably. To satisfy gauge invariance 
B and Bq should have the same normal component at the surface boundary. 

In the reduced MHD case we show that, even if it were possible to give a gauge 
invariant definition of the magnetic helicity, it would not have a physical meaning. In 
reduced MHD the field is decomposed as B = Bq + b^, where Bq = B^e^ is constant and 
uniform. We choose the vector potential Aq = B^xGy and b_|_ = V x A_|_ = V x (-i/^e^). 
We have only 4 terms which may contribute to an expression like (j2.93p . but only one of 
them is not null. In fact 

Ao • Bo = 5o a^e^^ • Bo e. = (2.94) 
Ao • b_|_ = i?o X • b_|_ = i?o X by (2.95) 
A_L • b_L = V' • b_L = (2.96) 
Ax • Bo = 5o ^ (2.97) 

Integrating over the volume equation (|2.95p vanishes because of the periodic boundary con- 
ditions of by, so only equation (j2.97p does not vanish. On the other hand, when integrating 
over the volume, this term is the zero frequency component of i.e. a constant which can 
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always be subtracted through a gauge transformation. 

In full ideal MHD magnetic helicity is a conserved quantity, but in reduced MHD 
it does not seem to have a physical meaning. While in 2D MHD magnetic helicity is zero 
but any moment of ip is conserved, in reduced MHD this does not apply. We rewrite for 
convenience equation ()2.46p with [(p,ilj] = ~V • (^u_|_): 

^._V.(,uJ+..| + J-Vi*. (2.98) 

When the z derivative term is zero, this equation becomes the 2D MHD equation, and it 
easily follows that any moment of is conserved when TZ^n — >■ oo. In the reduced MHD 
case the divergence gives a vanishing contribute (as it does for 2D MHD): 

/ d^x V • (V- u_l) = (ldaipu±-h (2.99) 
Jv Js 

because the normal component of the velocity is zero at the photospheric surfaces, and the 
remaining boundary surfaces are periodic. The presence of the z derivative term breaks the 
conservation of ip", except for a = 1, but in this case it is again a function which can be 
subtracted from ip through the gauge transformation 

^|J' = i^-VA fdt^ (2.100) 



so that for TZm — )• oo 



/ d^x V = const (2.101) 
Jv 



2.2 Anisotropic MHD Turbulence 

Wherever a fluid is set into motion turbulence tends to develop. When the fluid is 
electrically conducting, turbulent motions are accompanied by magnetic field fluctuations. 
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Although plasmas are abundant in the universe (it is said that 99% of the baryonic material 
in the universe is in the plasma state), conducting fluids are rare on earth, where electrical 
conductors are usually solid. Hence it is not surprising that magnetohydro dynamic turbu- 
lence (Biskamp [5j) has received attention only recently, after that hydrodynamic turbulence 
has been studied at length (Prisch [29]). 

A milestone in turbulence theory is the work by Kolmogorov [U] on the scaling 
properties of hydrodynamic turbulence, where he finds the well-known k~^^^ power spec- 
trum for total energy. The presence of the magnetic field strongly affects the properties of 
turbulence. While it is possible through a Galilean transformation to subtract a mean (or 
local) velocity field, this transformation has no effect at all on the magnetic field (a mean 
global or a local one). 

Since the pioneering work of Iroshnikov [l5] and Kraichnan [JH] (hereafter IK) 
there has been a lot of debate on which are the main properties of MHD turbulence. The 
Alfven effect, which arises from the fact that only oppositely propagating Alfven waves 
interact, and the hypothesis of homogeneity and isotropy lead to a k^^/'^ scaling for the 
energy spectrum, which differs from the Kolmogorov scaling in the hydrodynamic 

case (Kolmogorov [17], Obukhov [55]). 

The anisotropy of MHD turbulence is one the properties characterizing the recent 
debate (Shebalin et al. [84J, Sridhar & Goldreich [86], Goldreich & Sridhar [33j, Montgomery 
& Matthaeus [60], Goldreich & Sridhar [Mj, Cho & Vishniac [llj, Maron & Goldreich [5i] . 
Cho et al. [lO]). There is broad agreement that the anisotropy of MHD strongly affects 
its properties, simply due to the presence of a magnetic field, and that the hypothesis of 
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homogeneity and isotropy must be relaxed. 

Shebalin et al. [84] have shown that the energy cascade takes place mainly in the 
plane orthogonal to the static (DC) magnetic field, while it is weaker in the parallel direction. 
Sridhar & Goldreich [86] and Goldreich & Sridhar [33\ [M] have shown that the anisotropy 
gets stronger at large wave-numbers, i.e. whilst the cascade takes place. These results have 
been numerically investigated and confirmed by Cho & Vishniac [llj, Cho et al. [lOj . 

2.2.1 Turbulent Cascade and Phenomenology of the Inertial Range 

A characteristic property of fully developed turbulence is the presence of a broad 
range of different scales. Relevant physical quantities, such as energy, are excited within 
a certain spectral range k ~ called the injection scale. Nonlinear interactions transfer 
these quantities in A;-space towards larger wavenumbers up to the dissipation scale k ^ k^, 
where a dissipative physical process is supposed to act as a sink for this energy fiux. 

The region in Fourier space between the injection and the dissipation scale 

kin<.k<. kd (2.102) 

is called the 'inertial range'. In this spectral range the turbulence develops solely under the 
influence of the internal nonlinear dynamics without being directly influenced by either the 
external injection of energy or by dissipative processes. Spectra in the inertial range exhibit 
power-laws, and the inertial range can be defined as the wavenumber range within which 
the spectrum has a power-law behavior. 

What we have just described is called a direct cascade, but sometimes a fiux of 
energy in the opposite direction occurs, i.e. from the injection scale k ~ kin towards smaller 



40 

wavenumbers. This process is called inverse cascade, and when it occurs a second inertial 
range, besides (I2.102p . is found 

L-^ <^k<^ kin (2.103) 

The lower limit L^^ is usually determined by the size of the system L. 

We now briefly summarize the Kolmogorov phenomenology of the inertial range 
(K41). We explicitly consider the case of a direct cascade and of an isotropic system. The 
results that we will describe in the next chapter depart substantially from the K41 theory, 
but many concepts, ideas and notations are used also in anisotropic turbulence theory. 
Considering isotropic turbulence we can define the angle-integrated spectra in the following 
way: 

Ek= [ dn^ Ek, E = [ dk Ek (2.104) 



J Jo 

where thanks to isotropy we consider only the scalar wavenumber k = ^Jk^ + ky + k^. 

The dynamics of turbulence is controlled by the rate ej„ at which energy is injected 
into the system at the injection scale kin, it is subsequently scattered along the inertial range 
with the transfer rate e^, and finally swept away from the system at the dissipation scale 
k ^ kfi with the dissipation rate e^- For stationary turbulence all these spectral energy 
fluxes are equal 

em = et = ed = e. (2.105) 

This equality still holds approximately when the injection energy rate changes in time, 
because the more rapid dynamics at the small scales in the inertial and dissipation ranges 
adjust the spectrum rapidly compared to the slower dynamics of the large scales. For 
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convenience we divide the inertial range into a discrete number of scales In = 

lo>h> ••• > ^Af, i-e. ko < ki < ••• < /cat, (2.106) 

and the division is taken on a logarithmic scale /„ = 2~" Iq, where Iq is of the order of the 
large scale L. The time taken for the transfer of energy between two neighboring scales In 
and In+i is given by t„, the so-called eddy turnover time of the eddy dvi^, for simplicity 5vn 

T„ ~ ^ (2.107) 

dVn 

Since the energy flux e is constant across the inertial range, we can write 

e 1 —Sl (2.108) 

"^n '■n 

From this equality we can find the scaling 

6vn ~ 6^/3 (2.109) 

To obtain the energy spectrum we identify the eddy energy with the band-integrated Fourier 
spectrum 

rkn+i 

,2 



6vi^En^ / dkEk^knEk^ (2.110) 
from which we obtain, substituting kn A:, the —5/3 Kolmogorov spectrum 

Ek ^ e^l^ k-'^l^ (2.111) 



We remark again that the hypothesis of isotropy has been essential to obtain the 
K41 spectrum (|2.111l) . 
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2.2.2 The Alfven Effect 

In the hydrodynamic case isotropy is normally justified, but for plasmas a mag- 
netic field is always present. This introduces an anisotropy of the system which cannot be 
eliminated. 

Four decades after Iroshnikov |45j and Kraichnan [48j presented their ideas on 
MHD turbulence, the debate on which physical mechanisms drive it is still active. Rewriting 
the equations of incompressible MHD (j2.ip - (j2.5|) expressing magnetic field in velocity units, 
i.e. B — )• b = B/-^47rpo, and making them non dimensional choosing a characteristic velocity 
u* , a characteristic length I* , and the related crossing time t* = I* /u*, we have in terms of 
the Elsasser variables = u±b: 

^ = - (z^ • V) z± - VP + i V2z± (2.112) 
V-z=^ = (2.113) 

where we assumed the kinetic (TZ) and magnetic (TZm) Reynolds numbers, defined as 

1 v 1 n(? 

— = — > = > (2.114) 

TZ po ru* TZm 47r/Jo '*u* 

to be equal {TZm = TZ). P is the total pressure P = p/po + B^/Stt/jo in dimensionless form, 
and it is tied to the z^ fields by incompressibility (|2.113p : 

V^P = - ^ [d^zj) (2.115) 

In the following analysis we ignore the dissipative terms, because dissipation provided by 
TZ, supposed it has a high value, takes place only at small spatial scales. A constant and 
uniform magnetic field in absence of fluid motions: 

uo = 0, ho = VAe^ (2.116) 
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is a homogeneous solution of equations (j2.112p - (j2.113p . A linear analysis of incompressible 
MHD shows that we have only Alfven waves, in particular linearizing equations (12.1120 - 
(|2.113|) with the equilibrium (|2.116p yields: 

^TVA^ = 0, V-z± = 0. (2.117) 

Equations (|2.117p show that z~ describes Alfven waves propagating toward positive z at 
the speed Vj^, and z+ describes Alfven waves propagating at the same speed but in the 
opposite direction. A notable property of the Elsasser fields in equations (j2.112p - (|2.113p is 
the absence of self- coupling in the nonlinear term. In fact there is only cross- coupling of 
z+ and z . This property allows for a nonlinear generalization of the linear Alfven waves 
described by equations p.llTj) . Singling out the equilibrium (I2.116P (z^ = ±v_a in terms 
of Elsasser variables) 

z^=zf±VAe„ (2.118) 

and requiring the generalized Alfven wave z^ to retain its transversality, i.e. z^ -e^ = 0, we 
have for the nonlinear term 

9z^ 

z"^ -Vz^ = zf -Vzf ^VA^ (2.119) 

oz 

This means that if one of the two Elsasser fields zf is zero the nonlinear term in equa- 
tion (j2.119p vanishes and equations ()2.112p -( j2.113p assume the linear form (j2.117p for the 
remaining Elsasser field, even if its amplitude is not small. If at some time, like t = 0, 
z^ = and z^ = ¥{x,y,z) from equations ()2.117p we have that the solution is zf = 0, 
z^ = i{x,y, z — VA i)- If at time t = we had z^ = and zf = {{x, y, z), then the solution 
would be zj~ = and zf = f(x,y,2; + VAt). These nonlinear solutions are Alfven wave 
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packets of arbitrary form propagating nondispersively in the direction of the main field 
bo = vj[ez, and in the opposite direction. The dynamics are very simple as long as there is 
no spatial overlap ("collision") between two oppositely moving packets z+ and z . Hence 
only Alfven waves propagating in opposite directions along the guide field interact. This 
is the basis of the Alfven effect introduced independently by Iroshnikov [45j and Kraich- 
nan [48| who noted that the cascade of energy in MHD turbulence occurs as a result of 
collisions between oppositely directed Alfven wave packets. This result is quite general for 
MHD, in fact the guide field need not be an external static field, but can also be the field 
in the large-scale energy-containing eddies. 

2.2.3 The Iroshnikov-Kraichnan Formulation 

To derive the results of IK theory we consider a statistically steady, isotropic 
excitation of amplitude 5zf^ <C v_a, at the injection scale / of the equilibrium defined by 
equations (|2.116p . The turbulent cascade produces Alfven wave packets at scales A < / 
traveling in opposite directions along the large-scale field. A fundamental hypothesis made 
by both Iroshnikov |45j and Kraichnan |48j is that the energy transfer in Fourier space is 
local and isotropic. At this point we restrict the discussion to weak velocity-magnetic-field 
correlation 5z^ ~ 5z^ ~ Szx, which is the condition that applies to our model for a coronal 
loop and that we will discuss in more detail in the next chapter. 

We distinguish between two important dynamical time scales, the time for distor- 
tion of a wave packet 6z^ at scale A by a similar eddy Sz^, i.e. the eddy turnover time 

rx = T^, (2.120) 
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and the Alfven time t_4 = X/v^, which is the interaction time of the two oppositely moving 
wave packets. In general r_4 ^ r^, so the interaction time of the two wave packets is much 
shorter than the non- magnetic eddy turnover time tx- The change in amplitude A5zx is 
small during a single collision of two wave packets since it is proportional to the interaction 
time: 

~ — < 1. 2.121 

During successive collisions these perturbations add with random phases and, given the 
diffusive nature of the process, the number of collisions for the small perturbations to build 
up to order unity (i.e. A6zx ~ 5zx) is 

The energy-transfer time Tx, which in hydrodynamic turbulence is just tx, is longer 

Tx^ NxTA^^. (2.123) 
Making the substitution tx — )• Tx in equation (j2.108p we obtain for the spectral energy flux 

Ex ^4 T_A 



Tx A2 • 



E\ 5z\ TA 

e~-^ (2.124) 



From this we get the scaling 

5zxr^{evA\f'\ (2.125) 

and identifying the eddy energy with the band-integrated Fourier spectrum 5z\ ^ k Ek 
(where k ~ A^^) we obtain the Iroshnikov-Kraichnan (IK) spectrum for MHD turbulence 

Ek^{evA)'/'k-^/'', (2.126) 
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which is less steep than the k Kohnogorov spectrum (j2.11ip . Using the scahng (j2.125p 



in (I2.122P we obtain for the number of cohisions per energy transfer time: 




1/4 



(2.127) 



As we proceed along the cascade toward smaller scales A the number of collisions required 
for the fractional perturbations to build up to order unity increases, verifying the hypothe- 
sis (j2.12ip that we made at the beginning. Furthermore, during each collision the fraction 
of energy that cascades gets smaller with decreasing scale; in fact, from equations (|2.12ip 
and ()2.122p we have that 



In this sense we say that the cascade "weakens" at large wavenumbers. 

2.2.4 Beyond IK: Fully Anisotropic MHD Turbulence 

Isotropy is the underlying hypothesis used in the IK derivation of turbulence scal- 
ing properties (j2.120p - (j2.128p . In particular we have imposed the condition that the wave 
packets are isotropic on the length-scale A, also along the direction of the equilibrium mag- 
netic field bo = Vj^e^- 



Only later it has been understood that the anisotropy due to the presence of the 



main axial field not only acts through the Alfven effect, but has also a deep impact on 
the nonlinear dynamics, producing two different behaviors along the direction of the main 
field and in the orthogonal plane (Shebalin et al. [H], Sridhar & Goldreich [86], Goldreich & 
Sridhar [33], Montgomery & Matthaeus [60], Goldreich & Sridhar [34J, Cho & Vishniac [llj, 
Maron k Goldreich [M], Cho et al. [TO]). 



A5zx 



1 



(2.128) 
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Shebalin et al. [H] used the reduced MHD equations ()2.46p - ()2.47p to numerically 
investigate the cascade properties of a 2D turbulent system embedded in a strong field 
B = i?o ^a; directed along the x-axis, considering the system invariant along the z-direction, 
so to perform 2D numerical simulations in the x-y plane. As initial conditions they con- 
sidered wave packets with an isotropic spectral distribution in the k^-ky plane, as in the 
IK theory, so that the isocontours of the spectral densities were circles at the beginning 
of the simulation. They found that the spectrum evolves anisotropically by transferring 
energy to modes perpendicular to B far more rapidly than to modes with k parallel to B. 
In this way the initially circular spectral density contours elongated in the perpendicular 
direction. Even if it was not initially valid, the evolution proceeded toward the reduced 
MHD approximation. Even if the simulation was started with an isotropic initial condition, 
the temporal evolution was strongly anisotropic. 



The isotropic hypothesis used for the IK phenomenology (j2.120p - (|2.128p is there- 



fore neither consistent nor correct, as pointed out by Sridhar & Goldreich [86] and Goldreich 

& Sridhar j34j. In this sense Iroshnikov [l5] and Kraichnan [381 have only partially imple- 
mented the consequences of anisotropy in MHD turbulence through the Alfven effect, but 
the full consequences of anisotropy have been understood only later, and its elucidation is 
not yet complete. 

To understand why perpendicular transfer is easier we present a simplified pertur- 
bative argument in order to estimate possible energy transfer between Alfven wave modes. 
Introducing in incompressible MHD equations ()2.112p - (j2.113p the expansion 





(2.129) 
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where = itu^^e^ is the homogeneous equihbrium (j2.116p . and zf are Alfven waves of 
the form 

zf = ^A^exp[i(k-x±a;t)], (2.130) 

k 

where k • = 0, and uj = Vj^kz- Noting that from equation (I2.115P the pressure gradient 
term has only second order terms, i.e. 

VP ~ VPa + O (e^) , (2.131) 

at the second order we have 

This equation has basically the same structure of the wave equation for zf , except for an 
effective driving term on the right, due to the linear Alfven waves (j2.130p . 

The most efficient mechanism for fast energy transfer between modes is resonant 
interactions occurring among triads of modes with wavenumbers ki, ka and ka related by 
the conditions 

ki + k2 = ks, and cji + = (2.133) 

where ojj = f^k^^^. Shebalin et al. [HI] noted that the only nontrivial solution of (j2.133p 
requires that the z-component of one member of the triad, e.g. ks^^, must be zero. This 
implies that waves with values of not present initially cannot be created during collisions 
between oppositely propagating wave packets. Hence there is no parallel, i.e. along z, 
cascade of energy. Energy will cascade to large wavenumbers in the orthogonal plane k_L, 
making the turbulence cascade anisotropic. 
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We can now derive the anisotropic version of the IK theory (j2.120p - ()2.128p . taking 
into account that there is no cascade along the direction of the main magnetic field bo = 
va ■ We suppose again that the system is excited at the scale / in a statistically steady and 
isotropic fashion such that 6zi <^ v^- The absence of a parallel cascade implies that wave 
packets belonging to the inertial range have parallel scales / and perpendicular scales X < I. 
As previously supposed by IK, the Alfven effect takes place and only counter-propagating 
wave packets interact. The wave packets are "long-lived" and they need many collisions to 
loose a significant amount of energy. We distinguish again between the eddy turnover time 
characterizing the cascade in the orthogonal plane, 

TA = T^, (2.134) 
dzx 

where A ~ /cj^, and the Alfven time t_4, which is the interaction time of two oppositely 
propagating wave packets 6z^ and 6z^. Now, because of the absence of cascade along z, 
the Alfven time is scale-independent, i.e. it does not depend on the scale A: 

TA = — . (2.135) 

VA 

The interaction time is still small compared to the eddy turnover time r_4 ^ tx, so the 
energy loss of the eddy at the scale A is small during a single collision 

^ ~ ^ « 1 (2.136) 

Szx Tx 

The number of collision that a wave packet at the scale A must suffer for the fractional 
perturbation to build up to order unity is now 
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The energy-transfer time Tx is again given by 



while for the energy flux we have 



Tx ~ NxTA ~ (2.138) 



^2" — ■ (2.139) 

Tx VA 

From the previous equations we obtain the following scaling relation 

Szx-^[^f (2.140) 

and identifying the eddy energy with the band-integrated Fourier spectrum Szj^ ~ k_\_Ek^ 
(where ~ A) yields the anisotropic version of the IK spectrum for MHD turbulence 

E,^^{^^fkl\ (2.141) 

which exhibits the characteristic —2 spectral index. Another difference with the IK formu- 
lation is given by the behavior of the number of collisions at small scales: 

Nx-i'^Y X-i'-A^i (2.142) 



el^ J \6ziJ I 

Contrary to IK, Nx decreases with decreasing A. When the number of collisions decreases 
at small scales we say that the turbulence "strengthens". At a small enough scale the 
conditions (|2.136p - (|2.137p will not be satisfied, thus limiting the spectral range in which the 
spectrum Ej^^ ~ A;J^ applies. 

Weak perturbation theory (Zakharov et al. [96]) deals with the effects of the non- 
linear terms in equations (j2.112p - (j2.113p in a systematic, perturbative manner. When the 
nonlinear terms are ignored, the Fourier amplitudes and phases of the waves are constant 
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in time. However, the nonlinearity makes the amphtudes change slowly over many wave 
periods. It is this slow change in the amplitudes that determines energy transfer among 
the linear modes. A kinetic equation for the rate of change of energy in a mode with wave- 
vector k describes how other modes in the system affect the energy in this mode. To lowest 
order in the nonlinearity, the kinetic equation takes into account the interactions among 
modes taken three at a time, as shown in (I2.133p . When equations (12.133P is satisfied, one 
says that a 3- wave resonant interaction is allowed. When 3- wave resonant interactions are 
forbidden, because the 3-wave resonant coupling coefficients vanish, the effects of 4-wave 
resonant interactions must be considered. 

Goldreich & Sridhar [34J showed that although the anisotropic IK formulation (j2.134p - 
()2.142p describes a weak turbulence, in the sense that the fractional change in wave ampli- 
tude during each wave period is small, strains in the fluid are so strong that perturbation 
theory diverges. It turns out that not only 3-wave interactions contribute, but also higher 
order terms. At lowest order in perturbation theory, wave packets move along field lines. 
Thus the breakdown of perturbation theory can be understood physically by studying the 
geometry of the divergence of a bundle of field lines. Assume that the mean field lies along 
the z direction, and consider wave packets with longitudinal scale /, and transverse scale A, 
with A < /. For the turbulence to be weak we require the condition (j2.136p to be satisfied, 
i.e. 

X = ^ = ^ « 1 (2.143) 

so that N\ ^ 1. The rms inclination of the local field is 0\ ~ ^zxjvjs^^ so that after a single 
collision between two wave packets taking place along the longitudinal scale the wave 
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packet will suffer an orthogonal displacement S ~ I6x. Given the diffusive nature of the 
process after n collisions the wave packet has traveled a distance z ~ n Hn the longitudinal 
direction, suffering an rms orthogonal displacement 

A2 ~ n (5^ ~ nl'^el ~ \z\Wl. (2.144) 

The distance along z over which A increases by a factor of order A, i.e. ~ A^, is 

The perturbative expansion converges if the energy spectrum is cut off at small wavenum- 
bers, below A;^ L^, ~ 1 (Sridhar & Goldreich [86J, Goldreich &: Sridhar [S]). In this case it 
is shown that 3-wave resonant contributions vanish and 4-wave interactions must be con- 
sidered. For a system with a finite longitudinal extension L, such as a coronal loop, this 
means that for a sufficiently weak perturbation 6z\, we move from the anisotropic IK phe- 
nomenology (|2.134p - ()2.142p to a new one based on 4-wave resonant interactions that we 
now briefly describe. The elementary interactions involve scattering of two waves: 

ki -I- k2 = ks -I- k4, ui + UJ2 = u}3 + ooa. (2.146) 

Using ujj = vy[ kj^z and the z-component of the equation for k conservation, Sridhar &: 
Goldreich |86j proved that this scattering process leaves the fc^'s components unaltered. This 
implies that waves with values of that are not present in the external stirring cannot be 
created by resonant 4-wave interactions. The absence of a parallel cascade implies that wave 
packets belonging to the inertial range have parallel scales /, and perpendicular scales X < I. 
The phenomenology of this new cascade based on 4-wave interactions is very similar to the 
anisotropic IK phenomenology describe by equations (j2.134p - (j2.142p . and can be derived 
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in the same way with a few modifications. The Alfven effect takes still place, so that only 
counter-propagating wave packets interact, and the wave packets need many collisions to 
loose energy significantly. The eddy turnover time tx and the Alfven time r4 are defined in 
the same way (see equations (|2.134p and (|2.135p ) and have the same meaning, but now the 
fractional loss is different from what computed in (j2.136p . and from perturbation theory 
(Sridhar &: Goldreich ^86j, Goldreich &; Sridhar [34]) we have that: 



The number of collision that a wave packet at the scale A must suffer for the fractional 
perturbation to build up to order unity is now 




(2.147) 




(2.148) 



so the energy-transfer time Tx is 



Tx ~ NxT_A 



(2.149) 



3 ' 



and for the spectral energy fiux we obtain 




(2.150) 



Prom the previous equations the scaling relation 




(2.151) 



or equivalently 



2 




(2.152) 
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are obtained. Identifying the eddy energy with the band-integrated Fourier spectrum 
Szj^ ~ k±Ekj^ where A;_l ~ A we obtain the spectrum for MHD turbulence bases on 4-wave 
interactions: 

Ek^^e^^k'k (2.153) 
As in the case of the anisotropic IK phenomenology (|2.134p -( j2.142p the number of collisions 

decreases at small scales and the turbulence becomes stronger. 

The previous scalings (j2.134p - (|2.142p and (j2.147p - (|2.154p are based respectively on 
3- waves and 4- waves resonant interactions. These scalings can be generalized to the case of 
single collisions of n-waves. The derivation is very similar to what we have already done, 
and we briefly describe it. We suppose again that the system is weakly excited {5zi <^ vX) 
at the scale /, and that n- waves resonant interactions produce a perpendicular cascade. 
The interaction time is small compared to the eddy turnover time r4 <^ tx, so that the 
energy loss of the eddy at the scale A is small during a single collision, and is given by (see 
Goldreich & Sridhar |34j ) 

A6zx 



n 



5zx 

The number of collisions is now 



n-2 

~ ( ^ ) < 1, where n > 3. (2.155) 



/ r \ 2 / \ 2(n-2) / X \ 2{n-2) 



the energy-transfer time T; 

T^ ~ N^TA ~ 



Tx ~ NxTA - : ^;2„_5 , (2.157) 



-] ■ (2-158) 



6zx ~ e^I^ ( ^ j A— (2.159) 
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and the spectral energy flux 

From the previous equations we obtain the following scaling relation 

2n — 5 

or equivalently 

Tl-2 

and the anisotropic spectrum for MHD turbulence based on single n-waves resonant inter- 
actions is 

The spectral index spans from —2 and —7/3 for the cases that we have already treated 
(n = 3 and n = 4), and has a lower limit of —3 as n — )• oo. A common feature for the 
number of collisions 

N. . {£ (-) (2.162) 

is the "strengthening" of the turbulence at small scales, for all n > 3. 

A parameter that characterizes weak turbulence is x (see eq. (|2. 14311 ): 

X=T = — , (2.163) 

Aw^ Tx 

which, for weak turbulence, is small x ^ 1- At a sufficiently small scale along the cascade the 
turbulence will get strong enough so that x ~ 1) and just Nx ~ 1 collision (see eq. ()2.156p ) 
with another wave packet of comparable size will result in a fractional change in wave 
amplitude of order one, i.e. A5zx ~ 5zx from equations (j2.155|) - (|2.156|) . The same result is 
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obtained if the perturbation at the injection scale I is strong enough, i.e. 6zi ~ so from 
eq. (I2.163P we have x ~ 1- In this case wave packets lose their identity after they travel 
one wavelength along the field lines. Consequently the eddy turnover time and the Alfvenic 
time are the same, 

k^5zxr^ kiiVA, (2.164) 

where A ~ /cj^, which is called a "critical balance" (Goldreich & Sridhar [SH]). Consider 
an eddy of dimensions / and l± along the directions parallel and orthogonal to the mean 
magnetic field. Because of the turbulent transfer to smaller perpendicular scales, l± shrinks 
and the eddy becomes more elongated, leading to sheet-like structures limited only by 
dissipation. The spectral cascade takes place mainly in the orthogonal plane with constant 
energy flux e across the inertial range 

e ~ — ^ ~ const. (2.165) 
A 

Combining this relation with equation (j2.164p we obtain 

k. ~ — k]'^ ~ k]'^ ^-^/^, (2.166) 

where we define the scale ^ = v^je. This relation shows also that this anisotropy increases 
toward smaller scales, the ratio 

gets bigger for larger A;_|_ (smaller A). 

For the spectrum the original K41 phenomenology (|2.107p - ()2.11ip is valid, with a 
slight modification to account that the cascade takes place in the orthogonal plane. But 
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this analogy is only formal, because now we discard isotropy; in particular, the Alfven effect 
takes place and the cascade occurs mainly in the orthogonal plane. The eddy turnover time 
Tx and the energy transfer time now are the same 

TA ~ (2.168) 

and since the energy flux e is constant across the inertial range, we can write 

e^^^M, (2.169) 
Tx A 

yielding the scaling 

Szx^e^/^X^/^. (2.170) 
So from the band-integrated Fourier spectrum in the orthogonal plane (i.e. with A ~ kj^) 

Szl^k±Ek^ (2.171) 

we recover the —5/3 Kolmogorov spectrum 

~ 6^/3 (2.172) 

but now only in the orthogonal plane. 
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Chapter 3 



The Numerical Code 



A numerical code, written in Fortran 90 and parallelized with MPI, has been de- 



veloped to solve the non-dimensional reduced MHD equations ()2.46p - ()2.47p . that we rewrite 
here for convenience: 



The computational domain is a parallelepiped {i x i x L) with an orthogonal (x, y) square 
cross section of size ^, and an axial {z) length L, with the normalization 1 = 1 (because of our 
choice for the length-scale to render equations dimensionless shown in § I2.1.ip , and L > 1 
(see Figure [2?T]) . In the orthogonal planes [x and y directions) periodic boundary conditions 
are used coupled with a Fourier pseudo-spectral numerical method (Canuto et al. [8j). In 
the axial direction z, a velocity pattern is imposed at the top and bottom boundaries 
(see equations (j2.60p - (|2.6ip ). and a central finite difference scheme of the second order is 
used. Time is discretized with a third-order Runge-Kutta method coupled with an implicit 




(3.1) 



(3.2) 
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Crank-Nicholson scheme for the diffusive terms. 

When investigating turbulence, the diffusive terms provide a sink for the flux of 
energy at small scales. One of the problems present when performing a numerical inves- 
tigation is the limitation on the number of grid points. While the Reynolds numbers TZ 
and Tim have large values in most physical problems of interest, the number of points is 
necessarily limited when performing a numerical simulation which, in turns, restricts the 
Reynolds numbers to low values. Thus, diffusion is not limited to small scales, but also 
affects the large scales. If one of the purposes of the numerical investigation is to study the 
inertial range behavior, even with a resolution of 512 x 512 points in the planes, the inertial 
range is disturbed by diffusion. In numerical studies of turbulence it is often practical to 
use higher-order diffusion operators, or hyperdiffusion: 

^Vi ^ i_ Vi", a>l, (3.3) 

where the exponent a is called dissipativity. In this way dissipation is strongly concentrated 
at the small scales, and the dissipative pollution of the inertial range is avoided. We 
have performed simulations with both normal diffusion (a = 1) and hyperdiffusion with 
dissipativity a = 4. 

3.1 Fourier Transform and Spatial Derivatives 

Suppose h{x) is a complex periodic function of period £, defined on a uniformly 
spaced grid of N points: 

Xk = Ak, k = 0,...,N -I, where A = £/N. (3.4) 
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Indicating with /i^ the value of the function at the point x^, i.e. = h{Ak), the discrete 
Fourier transform i7„ at the wavenumber fc„ can be defined as (see Numerical Recipes |77j): 

N-l 

Hn =J2^k e-^^^"'^/^, n = 0, . . . , iV - 1, with fe„ = lirn/e, (3.5) 

k=0 

to which corresponds the inverse Fourier transform 

7V-1 

^'^ = ]^ J^i^ne*'""'/^, A; = 0,...,iV_i. (3.6) 

n=0 

Since the only differences between ()3.5p and (j3.6p are changing the sign in the exponential 
and dividing the result by N, a routine for calculating the discrete Fourier transform can 
also, with slight modifications, calculate the inverse transform. 

Although equation (13. 5p . the discrete Fourier transform, seems to be an 0{N'^) 
process (i.e. to compute the values of the function Hn it would require to compute A^ 
operations), an efficient algorithm called the fast Fourier transform or FFT [Numerical 
Recipes [77J) requires only C'(Alog2 A) operations. The FFT algorithm became generally 
known in the mid-1960s, from the work of Cooley and Tukey, and nowadays it is broadly 
used in scientific computing. FFTs are generally available as library subroutines and the fo- 
cus is mainly on achieving the best possible performance on a computing platform. We have 
chosen to use the FFTW library (Frigo & Johnson [2S|, see also http://www.fftw.org), 
which is Free Software distributed under the GNU General Public License. FFTW is typ- 
ically faster than other publicly-available FFT implementations and is competitive with 
vendor-tuned libraries that are tuned to work efficiently on specific CPUs. In this way we 
have achieved our goal to have a cross-platform portable code with a good performance. 

From equations (|3.4p and (j3.6p . it easily follows that the spatial derivative /i'^ 
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computed at the point j;^ is given by 

K = ^ (X.) = ^ E ^" ^ = 0' • • • ' ^ - 1- (3-7) 

n=0 ^ ^ 

Hence, to compute the derivative we first compute the Fourier coefficients Hn, then we 
multiply them by the factor i27rn/£, and finally we perform the inverse Fourier transform 
()3.6p with these modified coefficients. The precision of this pseudo- spectral method is much 
higher than that of an ordinary finite difference scheme. Even with relatively few grid 
points the computed value of the derivative is almost exact, with the error mainly due to 
the precision of the numerical processor (see Canuto et al. [8j). 

Time evolution is performed in Fourier space, rather than in coordinate space, so 
that we solve the Fourier transform of the reduced MHD equations (I2.46p - (l2.47p . Extending 
the notations introduced in ()3.4p . ()3.5p and ()3.6p to 2 dimensions, and noting that in our 
case £ = 1, we can write for a generic function /: 

f {x, y,z,t) = j^ Yl fr,s (z, t) e'C'r^+ksy) ^ ^here kr = 2^r, k^ = 2^s, (3.8) 

r,s=0 

and r, s are intergers ranging from to iV — 1. Introducing this expansion for the magnetic 
and velocity potentials, respectively and if, in ()2.46p - ()2.47p we obtain the reduced MHD 
equations in Fourier space: 

^ = ^ + [V^, i^l, - ^ K,s A,s, (3.9) 

^=^^^ + ^K,s^r,, (3.10) 

where k^ ^ = (27r)^ (r^+s^) and (r, s) / (0, 0). Note that to compute the Fourier components 
of the Poisson brackets, e.g. 

f dip dtp 
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we must first compute the inverse Fourier transform of (pr^s and V'r,s) then calculate the 
Poisson bracket in real space, and finally obtain its Fourier transform. 

When hyperdiffusion (|3.3p is used in Fourier space, the diffusive term is always 
linear, and the only thing that changes is the power of the A:-factor: 

-i^rs -J-^r"> ">!■ (3-12) 

Since Crank-Nicholson schemes are well-suited for the time advancement of linear terms, 
time evolution is performed with a third-order Runge-Kutta method (for the nonlinear terms 
and the z-derivatives) coupled with an implicit Crank-Nicholson scheme for the diffusive 
terms. To prevent numerical instabilities due to aliasing of the solutions, we truncate the 
Fourier transforms outside a circle of radius 2N/3 in the k^-ky plane, where is the number 
of modes and N x N the resolution in the real domain in the x-y plane (e.g. equation ()3.8p ). 
For time integration we require that the time-step satisfies the CFL (Courant-Friedrichs- 
Levy) condition in the x-y plane as well as in the z direction. This check is performed by 
a subroutine which dinamically adapts the value of the time-step. 

3.2 Message-Passing Interface and Parallel Computing 

To investigate some of the most interesting problems in physics, astrophysics, 
and engineering, challenging numerical computations are required. In numerical studies 
of turbulence we have already remarked that high grid resolutions are a necessity. For 
instance, if we had not performed some numerical simulations using hyperdiffusion ()3.3p and 
a numerical grid with 512 x 512 x 200 points {these are more than 52 millions grid points! 
52.428.800), most of the conclusions presented in this thesis would have been unattainable. 
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Serial computers are not suitable to perform this kind of simulation, and it is nec- 
essary to use a parallel computing system, i.e. a computer with more than one processor for 
parallel processing (commonly referred to as high-performance computers or supercomput- 
ers). These machines are used to perform numerical simulations of phenomena too complex 
to be reliably investigated by analytical methods, such as fully nonlinear problems, and/or 
very difficult or impossible to reproduce in a laboratory (unfortunately such as a coronal 
loop or most of the astrophysical problems of interest). 

There are many different kinds of parallel computers. One major way to classify 
parallel computers is based on their memory architectures. Shared memory parallel com- 
puters have multiple processors accessing all available memory as a global address space. 
They can be further divided into two main classes based on memory access times: Uni- 
form Memory Access (UMA), in which access times to all parts of memory are equal, or 
Non-Uniform Memory Access (NUMA), in which they are not. Distributed memory parallel 
computers also have multiple processors, but each of the processors can only access its own 
local memory. There is no global memory address space; but the processors communicate 
with each other through an intercommunication network, which can have many different 
topologies including star, ring, tree, hypercube, fat hypercube (a hypercube with more than 
one processor at a node), an n-dimensional mesh, etc.. Parallel computing systems with 
hundreds of processors are referred to as massively parallel. 

The Message Passing Interface (MPI) is a computer communication protocol (see 
[261 [271 [381 [2] , and also http://www.mpi-forum.org). The message-passing model posits 
a set of processes that have only local memory but are able to communicate with other 
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processes by sending and receiving messages. It is a defining feature of the message-passing 
model that data transfer from the local memory of one process to the local memory of 
another requires operations to be performed by both processes. Although the specific com- 
munication network is not part of the computational model it is one of the more delicate 
points in parallel computing, often representing the bottleneck to the performance of a nu- 
merical code. The challenge is always to try building intercommunication networks (also 
called switches) that keep up with speeds of advanced single processors. Faster computers 
require faster switches to enable most applications to take advantage of them. 

MPI is an attempt to collect the best features of many message-passing systems 
that have been developed over the years and to improve and standardize them. MPI is a 
library, not a language. It specifies the names, calling sequences, and results of subroutines 
to be called from Fortran, C and C++ programs. Programs are compiled with ordinary 
compilers and linked with the MPI library. It is emerging as a standard for communication 
among the nodes running a parallel program on a distributed memory system, although MPI 
can also be used on shared memory computers. Its advantage over older message passing 
libraries is that it is both portable (because MPI has been implemented for almost every 
distributed memory architecture) and fast (because each implementation is optimized for 
the hardware on which it runs). 

The message-passing model fits well on separate processors connected by a com- 
munication network. Thus, it matches the hardware of most of today's parallel supercom- 
puters. Where the machine supplies extra hardware to support a shared-memory model, 
the message-passing model can take advantage of this hardware to speed up the rate of 
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data transfer. Message passing has been found to be a useful and complete model in which 
to express parallel algorithms. It provides control when dealing with data locality and, by 
controlling memory references more explicitly than any of the other models, the message- 
passing model makes it easier to locate erroneous memory reads and writes. 

But the most compelling reason why message passing will remain a permanent 
part of the parallel computing environment is performance. As modern CPUs have become 
faster, management of their caches (divided in many levels) and the memory hierarchy 
in general has become the key to getting good performance. Message passing provides a 
way for the programmer to explicitly associate specific data with processes and thus allow 
the compiler and cache-management hardware to function fully. Indeed, one advantage 
distributed-memory computers have over even the largest single-processor machines is that 
they typically provide more memory and more cache. Memory-bound applications can 
exhibit superlinear speedups when ported to such machines and, even on shared-memory 
computers, the message-passing model can improve performance by providing more pro- 
grammer control of data locality in the memory hierarchy. 

For these reasons message-passing has emerged as one of the more widely used 
paradigms for implementing parallel algorithms. Although it has shortcomings, message- 
passing comes closer than any other paradigm to being a standard approach for the im- 
plementation of parallel applications. Message-passing has only recently, however, become 
a standard for portability. Before MPI, there were many competing variations on the 
message-passing theme, and programs could only be ported from one system to another 
with difficulty. 
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3.3 Code Parallelization 

The numerical code has been developed to solve equations (|3.9p - (j3.10p . which are 
the reduced MHD equations in Fourier space. As already noted in § 13.11 at each time- 
step the Fourier transforms of the fields and their inverse transforms must be computed to 
calculate the Poisson brackets. As shown in ()3.5p and (|3.6p to compute the value of a Fourier 
transform at one point along a given direction, we must know the values of the function in 
all grid points along that direction. This makes the Fourier transform an intrisically non- 
parallel computation. In fact if the values of the function along that direction were assigned 
to different processors, at each time-step all these processors would have to communicate 
these values to one processor which would compute the Fourier transform. Communications 
between processors is not very fast and then we want to minimize it. 

Although there are algorithms which parallelize the fast Fourier transform more or 
less efficiently, as a first step we have avoided this, taking advantage of the fact that we use 
a finite difference scheme of the second order along the z direction. These schemes are very 
well suited for parallelization, because to compute the value of the derivative of a function 
in one point we only need to know the values of the function in a few neighboring points. 
In the case of our central finite difference scheme of the second order only the values of the 
two neighboring points are required. We have decomposed our computational box, which is 
a parallelepiped (see Figure 12. ip of dimensions 1 x 1 x L with a grid of n x n x points, 
into slices along the z direction. So the grid points lying on a x-y plane at z = const are 
assigned to the same processor and no communication is required to perform the numerical 
computations in the plane, including the FFT. To compute the z-derivatives, where the 
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two neighboring points belong to the same shce, no communication is needed. This always 
happens except for the points at the top and bottom boundary x-y planes of a single 
slice. From each slice the values of the functions in these two boundary planes must be 
communicated, at each time step, respectively and exclusively to the processors to which 
have been assigned the next and previous slices of the computational box. 

The next task is to decide how to assign processes to each part of the decomposed 
domain. Handling this assignment of processes to regions is one of the services MPI provides 
to the programmer because the best (or even a good) choice of decomposition depends on 
the details of the underlying hardware. As noted in § 13.21 the processors of a supercomputer 
are linked by an interconnection network through which communication takes place. These 
networks can have many different topologies. Most of these topologies are complicated, 
and especially on massively parallel computers, when the code uses many processors it is 
very unlikely that each of these processors can communicate directly with all the others. 
Given our choice for the decomposition of the computational box, ideally we would like 
the processors which have been assigned to two neighboring slices of the computational 
box could communicate directly between themselves. In general the description of how the 
processes in a parallel computer are connected to one another is often called the topology of 
the computer (or more precisely, of the interconnection network). In most parallel programs, 
each process communicates with only a few other processes; the pattern of communication 
is called an application topology or virtual topology. 

It might seem that simply assigning processes in increasing rank from the bottom 
is the best approach. On some parallel computers, however, this ordering can degrade 
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performance. It is hard for anyone but the vendor to know the best way for application 
topologies to be fitted onto the physical topology of the parallel machine. MPI allows 
the vendor to help optimize this aspect of the programming through implementation of 
the MPI topology functions. MPI allows the user to define a particular application, or 
virtual, topology. An important virtual topology, the Cartesian topology, is a decomposition 
in the natural coordinate (e.g. z) direction. Although topology functions are sometimes 
treated as an exotic and advanced feature of MPI, they make many types of MPI programs 
easier to write. We implemented a one dimensional Cartesian Topology, and found it very 
useful. In this way MPI assigns the slices of the decomposed computational box to processes 
such that the topology of the underlying interconnection network and the topology of the 
communication in the numerical code match in the best possible way. 



3.4 Performance Evaluation 

Two quantities which are commonly used to quantify the performance of a numer- 
ical code on a parallel machine are the so-called speedup and efficiency. 

A parallel system is defined as the implementation of a parallel algorithm on a 
specific parallel computer. The dimension of the problem W is the number of operations 
required by the fastest known serial algorithm, and is equivalent to the concept of compu- 
tational complexity. 

We distinguish between the serial execution time Tg, the time interval between 
the beginning and the end of the execution of the program on one single processor, and the 
parallel execution time Tp, the time interval between the beginning of the execution and the 
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instant of time in which the last processor ends the execution. The serial execution time is 
basically a function of the computational complexity W, while the parallel execution time 
depends on W, the number of processors used Np, and the kind of interconnection network 
present in the supercomputer. The speedup S for a parallel system is defined as the ratio 
of the serial to the parallel execution time: 



It could be naively thought that if the execution time on a single processor is Tg , then the 
execution time on Np processors should be Tp = Tg/Np, in which case for the speedup we 
would obtain the linear behavior 



Unfortunately this relation is rarely found, and when it happens, the reason is not because 
Np processors perfectly distribute among themselves the work-load. 



overhead. The first, obvious, source is that when a numerical code is parallelized, it is 
usually necessary to change the structure of the code itself, adding some extra computations 
with respect to the serial version of the code. Other important sources of overhead are the 
communications and the imbalance of the work-load among the processors. But, especially 
for computations involving large data-sets (e.g. high-resolution simulations with a lot of 
grid points), parallelization has the advantage of splitting the computational box among 
many processes so that each process is assigned a relatively small number of grid points. In 
this way each single processor has to execute fewer computations, but - more importantly 




(3.13) 



S{W,Np) = Np. 



(3.14) 



In fact, as already partially remarked, there may be many different sources of 
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- since the data occupy a smaller amount of memory, the processors can manage their 
caches in a more efficient way. Sometimes the cache-management improves dramatically, 
leading to so-called superlinear speedups, i.e. to speedup values which exceed the linear 
behavior S {W, Np) > Np. 

In general the competition between the improved cache-management and the over- 
head due to the communication over the interconnection network determines the perfor- 
mance of a numerical code on a parallel system. 

EfHciency is defined as the ratio of the speedup (j3.13p to the number of processors 



Its value, typically included between zero and one (except for superlinear algorithms for 
which E > 1), estimates how efficiently the processors are used, compared to the time 
wasted in communication and synchronization. Algorithms with linear speedup (|3.14p and 
those running on a single processor have E = 1, while many algorithms that are difficult to 
parallelize have E ~ 1/logA'p that approaches zero as the number of processors increases. 

Figure [3TT] shows two speedup curves for our numerical code relatively to the SGI 
Altix 3000 supercomputer at Naval Research Laboratory. This computer has a total of 
64 Intel Itanium2 processors with 1.3 GHz of clock speed and 3 MB of cache memory per 
processor. They are interconnected with a Myrinet network. In most of the simulations we 
have performed we have typically used two different numerical resolutions, a lower resolution 
grid with UxXnyXUz = 256 x 256 x 100 points and a higher resolution one with Ux xuyXn^ = 
512 X 512 X 200 points. The speedup for the lower resolution case is shown on the left in 



p 




p 



(3.15) 
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Figure 3.1: Speedups of our numerical code on the SGI Altix 3000 computer at Naval Re- 
search Laboratory (NRL, Washington, DC). The computer has a total of 64 Intel Itanium2 
processors at 1.3 GHz of clock speed and 3 MB of cache memory per processor. The dashed 
line represents the linear behavior ()3.14p . On the left it is shown the speedup for simulations 
with a grid of n^. x Uy x = 256 x 256 x 100 points. On the right we show the speedup 
for simulations with a higher resolution grid with Ux x Uy x Uz = 512 x 512 x 200 points. 
The scalability of the numerical code increases at higher resolution. 



Figure [3Tl An almost linear behavior is shown up to twelve processors, while for an higher 
number of processors the performance starts to decrease. On the other hand using the 
higher resolution grid, the ratio between computations and communication grows favorably 
for computations, and the corresponding speedup curve, shown on the right in Figure [3T| 
actually exhibits an approximate linear behavior up to thirty processors. 

The property of a good performance over a wide range of processors is called 
scalability. We have achieved our goal of developing a cross-platform portable parallel code 
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with a good scalability. 
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Chapter 4 



Linear Analysis 



In this chapter we present the results of our numerical and analytical investigations. 



We make a Cartesian model (see §[2]) of a coronal loop, i.e. the loop is "straightened out" so 
that the computational box is a parallelepiped (see Figure 12. ip of unit square cross section 
and axial length L > 1 (with overall dimensions 1 x 1 x L). The dynamics of the plasma 
are studied using the dimensionless equations of reduced MHD (|2.46p - (|2.47p that we repeat 
here for convenience: 



where is the ratio between the axial Alfvenic velocity (associated to the strong, uniform, 
homogeneous field Bq = Bq Bz) and the typical rms photospheric velocity Uph 



94' dip . 1 , 



(4.1) 




(4.2) 



Bo 1 



(4.3) 



VA = 
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From observations (see § [I]) we know that the typical photospheric magnetic field is of the 
order of Bq ~ 10 gauss, the typical numerical electron density ~ 10^*^ cm^'^ (i.e. a mass 
density pq ~ 1.7 • 10^^^ gcm~^), and for the photospheric motions Uph ~ lkms~^, so the 
Alfven velocity is ~ 10^ kms^^, and the dimensionless ratio (j4.3p ~ 10"^. 

In (|4.1|) - (|4.2p we have supposed that the magnetic and kinetic Reynolds numbers 
are equal (TZm = T^)- The velocity and magnetic scalar potentials, respectively if and ip, 
are linked to the orthogonal components of the physical fields by 

ujL = V X {(fe^) , b_L = Vx(^e^), (4.4) 

and to the axial component of the electric current (j) and of the vorticity (u) by 

w = (V X u^)^ = -Vl^, j = (V X bx), = -ViV. (4.5) 

In the orthogonal x-y planes periodic boundary conditions are imposed, while in the z di- 
rection at the top {z = L) and bottom [z = 0) planes — that represent the two photospheric 
cross-sections to which a coronal loop is anchored — a velocity pattern is imposed, i.e. the 
velocity potential is specified. 

As initial conditions we specify the velocity patterns on the bottom and top 
planes ip{x,y,d) and ip{x,y,L), and v^. The potentials 99 and ■0 inside the computational 
domain are set to a small perturbation. The parallel numerical code (see § [3]) solves the 
Fourier transform of reduced MHD equations (|3.9|) - (|3.1U|) using normal diffusion or hyper- 
diffusion (see equations (j3.3|) and (j3.12p ) with dissipativity a = 4. 

In § S] we present the linear analysis of our system, which is characterized by the 
propagation of Alfven waves, while in § [H] the phenomenology of the turbulent dynamics 
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which develops in the nonhnear stage is discussed. Finally in §[5] the numerical simulations 
are presented. 



4.1 Boundary Conditions and Linearity 

The velocity forcing imposed at the boundaries produces Alfven waves that prop- 
agate into the computational box. The scalar fields ip and ip associated with these waves 
are small and since the Poisson brackets are quadratic in these terms — so their contribu- 
tion is negligible respect to the linear terms — the initial stage of the dynamics is linear. 
Linearizing the equations (|4.ip - ()4.2p . and neglecting for the moment the dissipative terms, 
yields 

d'^ dip dip dilj 

Introducing z^, the analogs of the Elsasser variables (j2.52p for the scalar potentials: 

z^ = ip±'4) (4.7) 



equations (j4.6p can be written as 



dz'^ dz^ 



at =^"-&r 

These wave equations show that z^ describes a wave propagating towards negative z, while 
z" waves propagating towards positive z. In terms of these Elsasser variables, the boundary 
conditions need be expressed only for the incoming wave, i.e. z~ at the bottom {z = 0) and 
z^ at the top {z = L). The imposition of the velocity potential patterns ip^ and p^ 
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respectively in z = and z = L is then achieved by the reflection condition: 



z \z=o = {-z^ + 2ip) |2=o = -^"'"12=0 + 2(/3°(x, y, t), 
z~^\z=L = {-z~ + 2ip) \z=L = -z~\z=L + 2ip^{x,y,t), 



(4.9) 
(4.10) 



Introducing the column-vector U and the square matrix A 



U 



z^ 



VJ 



A = VA 



I. 0^ 



equations (|4.8p can be written as 



dt dz 



indicating with C the operator: 



C 



— - A— 

dt dz 



(4.11) 



(4.12) 



(4.13) 



linear equations (j4.12p . or equivalently (14. Sh can be written as C{\J) = 0. Because of 
linearity, if Ui and U2 are solutions, then aUi + /3U2 is also a solution. If Ui, U2 have 
respectively the stream- functions ip\, ipi and ip2, (/^f boundary velocity potentials (j4.9p - 
(j4.10p . then a simple calculations shows that the boundary conditions satisfied by the new 
solution 



U 



are, indicating with ip^, ip^ the boundary stream- functions of the solution U: 



z 



aUi + /3U2 



az^ + (3z2 
^az^ + /3z2 J 



(4.14) 



ip^ = - [z^ + z 



[z-^ ^z- 



=0 = - [a (4 + zf) + /3 (4 + Z2-)]^^o = + fip>\ (4.15) 
=L = ^ [« (4 + ^r) + /? (4 + ^2")] ,=L = "V^f + ^H^k (4.16) 
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which is simply the hnear combination of the boundary conditions for the solutions Ui and 
U2. We will use the property (I4.15p - (l4.16p many times in the following paragraphs. In 
particular, if Ui is the bottom plane forced solution with cp^ = ip^ , ip^ = 0, and U2 is the 
top plane forced solution with = , y:i2 = then the function U = Ui + U2 is the 
solution that satisfies the boundary conditions 

(z+ + z-)|,=o = 2/, {z+ + Z-) U=L = 2ip^ . (4.17) 

This result allows us to investigate the two one-sided solutions separately, the two-sided 
solution results simply from their sum. 



4.2 One-sided Problem: Time Independent Forcing 

Considering at first only time-independent velocity patterns, calling r_4 = L/vj[ 
the crossing time of an Alfven wave for the axial length L, the solution to the linear equa- 
tions (j4.8p with the top boundary forcing 

{z+ + z-)\,=o = 0, {z+ + z-)U=L = 2ip^{x,y), (4.18) 

is given by 

for 2k Tj, < t < {2k + 1) fc = 0,1,2,... (4.19) 

{x,y,z,t) =2kip^ + 2ip^ (4.20) 

z- {x,y,z,t) = -2kip^ (4.21) 
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while 



for {2k + 1) TA<t < {2k + 2) TA k = 0,l,2,... (4.22) 

z+ {x,y,z,t) = 2kip^ + 2ip^ (4.23) 

{x,y,z,t) = -2kip^ -2ip^ (4.24) 

0^ and are step functions "propagating" towards negative and positive z, starting 
from one side and reaching the other one distant L in the time r4. For example, indicating 
with Q the step function 

1 Z > Zq 



(4.25) 

Z < Zq 



e{z-zo) = < 

referring to the time interval ()4.19p 

e^{z,zo) = &{z- zo) where zq (t) = L - VA^t - 2kTjs,) (4.26) 

so that the point zq moves from zq = L at time t = 2A; to = at time t = {2k + 1) r_4, 
at the speed va- In the same fashion, referring to the time interval ()4.22p 

6^ (z, zq) = G (zo - z) where zq {t) =VA[t- {2k + 1) ta] (4.27) 

so that the point zq moves from = at time t = {2k + 1)ta to zq = L at time t = 
{2k + 2) TA, at the speed va- 

To understand why the solution to the problem (j4.8p with the one-sided boundary 
condition (j4.18p has the structure shown in ()4.19p - ()4.24p . consider the first interval of time 
in equation (|4.19p < t < r^. At time t = the boundary condition z+ = 2(p^ in z = L 
causes a z+ wave to propagate towards negative z, while the boundary condition z~ = 
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in z = does not excite any z~ wave to propagate towards positive z. Hence the solution 
during this first time interval is: 



= 299^6^ , =0 (4.28) 

which are exactly equations (l4.20[ )-( l4.2ip for = 0. When at time t = tj^ the z'^ wave 
reaches the z = boundary conditions (|4.18p read: 

z-U=o = -2+U=o = -29?^, z+U=L = -z-\^=L + 2(^^ = +2^9^ (4.29) 



which cause the z"*" wave to continue propagating towards negative z, and a new wave 
to originate from z = (at time t = tX) and propagate towards positive z (reaching z = L 
at time t = 2t_a). This is described by equations (|4.23p - (j4.24p for k = 0, i.e. 

z+ {x,y,z,t) = 2^^ , z- {x,y,z,t) = -2ip^Q^ (4.30) 

The iteration of these arguments finally leads to equations (|4.19p - (j4.24l) . 

Solution (j4.19p - ()4.24p shows that "reflection" boundary conditions (j4.18p result in 
a wave propagation in our computational box that causes the amplitude of the scalar fields 
z^ to grow of the value 2ip^ at each interval of time corresponding to an Alfvenic crossing 
time r4. In terms of the velocity and magnetic potentials equations (|4.19p - (j4.24p yield: 

for 2k Tj, < t < {2k + 1) fc = 0,1,2,... (4.31) 

^{x,y,z,t) =^^e^ (4.32) 
ij{x,y,z,t) = 2kip^ + ip^e^ (4.33) 
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for {2k + 1) Tjs, < t < {2k + 2) Tjs, /c = 0, 1, 2, . . . (4.34) 

ip{x,y,z,t)=^^ -ip^Q^ (4.35) 
ij{x,y,z,t) = 2kip^ + ip^ + ip^e^ (4.36) 

From these equations we note that while the velocity potential ip is a. step function that 
bounces back and forth with the amplitude ip^, the magnetic potential is characterized 
by a propagating wave and grows roughly linearly in time gaining a factor ip^ at each 
time-interval r^. From equations (14. 4p we obtain for the velocity and magnetic fields: 

for 2kT_A<t< {2k + 1) A: = 0, 1, 2, . . . (4.37) 

ux(x,2/,z,t) =u^(a:,2/) (4.38) 

(x, y, z, t) = 2k (x, y) + (x, y) 9^ (4.39) 

and 

for {2k + 1) TA<t<{2k + 2) T_A A: = 0, 1, 2, . . . (4.40) 

(x, y, z, t) = (x, y) - (x, y) (4.41) 

(x, y, z, t) = {2k + 1) (x, y) + (x, y) 9^ (4.42) 

where (x, y) is the velocity forcing imposed at the top boundary plane z = L. Equa- 
tions ()4.37p - ()4.42p show that this velocity forcing results in a mapping of the boundary 
velocity inside the computational box for both the magnetic and velocity fields. The differ- 
ence is that while the velocity is bounded to its boundary value , the magnetic field does 
not satisfy this condition and grows linearly in time, its amplitude gaining a factor at 
each Alfvenic crossing time t_4. 
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4.3 Physical Interpretation 

The physical interpretation of equations ()4.37p - (|4.42p is simple. The photospheric 
motions imposed at the boundary plane cause the footpoints of the magnetic field-lines 
(which are almost straight because of the strength of the dominant axial magnetic field) to 
move. As a result Alfven waves, of amplitude equal to the velocity forcing u^, are launched 
into the computational box and reflected at the boundaries. Because of the boundary 
conditions the magnetic field grows roughly linearly in time, while the velocity field is 
bounded to its boundary value. Like in the Parker picture (see § \1.3.^) . the shuffling of 
the field-lines footpoints causes the field-lines to braid. But differently from that picture we 
make explicit from equations |^.371 )-( f^.^g[ ) that the braiding of the magnetic field-lines takes 



place through Alfven wave propagation. This is a phenomenon that naturally happens for 
a system embedded in a strong axial field, in both the linear and nonlinear stages. On the 
other hand this simple observation leads to one of the main original results of this thesis: 
since the field-line braiding is due to propagating Alfven waves then the nonlinear stage 
dynamics (which will be analized m § is described by anisotropic turbulence (see § 2.2.4\ ) 
due to the "collisions" of counter-propagating Alfven waves packets. 

4.4 Two-sided Problem: Time Independent Forcing 

In § 14.21 we have obtained the solution for the linearized equations (j4.8p with a 
velocity forcing at the top boundary 

(z+ + z-) |,=o = 0, {z+ + z-)\,=L = 2vL^{x,y), (4.43) 
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Using the same techniques, the solution with a velocity forcing at the bottom boundary 

(z+ + Z-) |,=o = 2u0, (z+ + Z-) l=L = 0, (4.44) 



is the following: 



for 2kTA<t< {2k + 1) A; = 0, 1, 2, . . . (4.45) 

ux(x,y,z,t) = u°(x,y) (4.46) 
(x, y, t) = -2k uO (x, y) - (x, y) e_, (4.47) 



and 



for {2k + 1) TA<t <{2k + 2) TA A; = 0, 1, 2, . . . (4.48) 

(x, y, t) = u° (x, y) - u° (x, y) 6^ (4.49) 

(x, y, i) = - (2A; + 1) u° (x, y) - (x, y) 6^ (4.50) 

And finally, the general two-sided solution with velocity patterns u° at the bottom boundary 
plane z = Q and at the top boundary plane z = L, is given by 

for 2kTA<t< {2k + 1)ta A; = 0, 1, 2, . . . (4.51) 

(x, y, z, t) = u° (x, y) 8^ + (x, y) 9^ (4.52) 

bx (x, y, z, = 2k [u^ (x, y) - u° (x, y) ] - u° (x, y) 0^ + (x, y) (4.53) 

and 

for (2A; + 1) < t < (2A; + 2) A; = 0, 1, 2, . . . (4.54) 

(x, y, t) = u° (x, y) + (x, y) - u° (x, y) - (x, y) 8^ (4.55) 

bx (x, y, t) = (2A; + 1) [u^ (x, y) - u° (x, y)] - (x, y) 0^ + (x, y) 0^ (4.56) 
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4.5 One-sided Problem: Time Dependent Forcing 

Solutions (j4.19|) - (|4.24|) can be generalized to the case of a time-dependent forcing: 

for 2kTA<t< {2k + 1) A: = 0, 1, 2, . . . (4.57) 

fc-i 



and 



/ L-z\ 
(x,y,z,t) = +^2(^^ it-2nTA ) 

n=0 ^ / 



+2 if^ (^t - 2kTj^ - (4.58) 
(x, y, z, t) = - V 2 (^^ f t - (2n - 1) - — ) (4.59) 



for {2k + l)TA<t < (2k + 2) T_A A; = 0,1,2,... (4.60) 
z+ (x, y, z, t) = + V 2 f t - 2nT^ - (4.61) 

z- (x, y, z, t) = - V 2 ( t - (2n - 1) - — ) 



-2 - {2k + 1)ta- — J (4.62) 

These equations describe the propagation and reflection of the Alfven waves originating 
from the boundary z = L. It can be easily checked that when the stream- function (p^ is 
time-independent equations (|4.57|) - (|4.62|) coincide with (|4.19|) - (j4.24|) . 

Most of the energy of photospheric motions is at low frequency (~ 3.3 mHz) and 
at spatial scales of ~ 1000 km. But a small fraction of their energy is also present at smaller 
spatial scales and higher frequencies. It is therefore interesting to perform an analysis of a 
Fourier component 

{x,y,t) = f {x,y)cos{ujt) . (4.63) 
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Substituting this expression in the previous equations (|4.57p - (|4.62p we note that, in general, 
the waves which are continuously injected and reflected in the computational box will be 
out of phase, so their sum will be of the order of /. The only exception to this behaviour is 
when they are in phase. This happens when the terms that are integer multiples of cj • 2 r4, 
arising in the cosine terms from the substitution (j4.63p in (j4.57p - (j4.62p and which are the 
only ones to differentiate among them the terms in the summations, are equal to an integer 
multiple of 2tt, so that they can be removed. The term uj ■ 2nr^ with n integer to be a 
integer multiple of 2tt requires the "resonant" condition 



to be satisfied. Using vji = 1/r^, the corresponding frequencies v = oj/Itt are given by 



(^m TA = rriTT 



with 



m G N. 



(4.64) 



■m 



m 1 



with 



m G N. 



(4.65) 



At the frequency cJ; 



rmr/TA equations ()4.57p - (|4.62p yield: 



for 2kT_A <t < {2k + 1) TA 



fc = 0,l,2,... 



(4.66) 




(4.67) 



(4.68) 



while 



for {2k +l)TA<t<{2k + 2) ta 



A; = 0,1,2,... 



(4.69) 




(4.70) 



(4.71) 
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In terms of potentials, remembering that cos (a it /3) = cos a cos /3 =F sin a sin /3 we have 



for 2kTj^<t <{2k^\) A; = 0,1,2,... 



LP (x, y, z, t) = -2 (-!)"" /c / (x, y) sin {urn t) sin w. 



+ (-1) / (X, y) cos ( t + 



tp (x, y, z, t) = +2 (-l)"" /c / (x, y) cos (w^ t) cos ujm — 

VA 



+ (-l)™ / (X, y) ©4- cos \ujrnt + ^^' 



VA 



(4.72) 



(4.73) 



(4.74) 



and 



for {2k + I) Tji^ < t < {2k + 2) Tj, A: = 0,1,2,... 



if (x, y, z, t) = -2 (-1)" A; / (x, y) sin (w^ t) sin a;„ — 



+ (-l)™/(x,y) 



cos \UJrnt + ^m 



VA 



0^ COS f UJm t — UJm 



VA 



■0 (x, y, z, t) = +2 (-1)™ A: / (x, y) cos (w^ t) cos ( — ) 



+ (-l)™/(x,y) 



cos [ (jJmt + U>m — ] + COS ( t - — 

vaJ V 



(4.75) 



(4.76) 



(4.77) 



The main result of these calculations is that, when present, these so-called resonances grow 
linearly in time with their amplitude proportional to time and to the forcing function ()4.63p 
f {x,y), i.e. to the amplitude of the boundary motions at the corresponding resonant fre- 
quency. This implies that as long as the energy fraction of the high frequency components is 
small, their contribution to the linear dynamics, and hence also to the nonlinear dynamics, 
is and remain small. 
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4.6 Effects of Diffusion 

In all the previous calculations we have neglected the effect of diffusion. In general 
diffusion is not very important when the Reynolds numbers are very high but numerically 
we use lower values so that it can affect the solutions that we have found (even at large 
scales). In particular diffusion is more important when the linear stage lasts for long times. 
The mathematical algebra is simplified if we consider the physical fields averaged in time, on 
scales bigger than the Alfvenic crossing time r^. For example from equations (I4.37p -( l4.42p 
we have 

u^'^u^{x,y) y, b^r^u^{x,y) — (4.78) 

While diffusion will only slightly change the shape in z of the velocity field in (j4.78p . it has 
a stronger effect on the magnetic field which would otherwise grow linearly in time. Here 
we describe briefiy such effect. Linearizing equation (j4.ip retaining the diffusive term, in 
terms of the magnetic and velocity fields we have 

db± du± 1 „9 , / , „ N 

We use the general boundary conditions u_l = u'^, = respectively in z = and 
z = L, but now we take into account that our forcing velocity has only components at 
the injection scale. In dimensionless form, the cross-section or our numerical box is a 
unit square of lenght i = 1, and the injection scale is im ~ 1/4. For example, one of 
the forcing patterns we use is = ey sin{27rr x + 1) with r = 4, corresponding to the 
wavelength A = = 1/4. For this pattern the relation V5_'U^ = — (27rr)^ is exact, 
but it is still approximately valid when we consider velocity patterns which result from 
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the linear combination of components with A ~ 1/4, i.e. r ~ 4. For the integration of 
equation (I4.79P we use the ansatz V^fo^ = — (27rr)'^ This is justified because the 
magnetic field is always the result of the mapping of the boundary velocity forcing, as 
found in the calculations without diffusion. Then integrating (j4.79p over z, dividing by the 
length L, using the boundary conditions u_|_(z = 0) = u'^ and u_|_(z = L) = u^, we obtain 
for the field 6^ averaged in z: 



(4.80) 



and indicating with Up/j = — u'^, and with t-ji = TZ/ (27rr)^ the diffusive time-scale and 
with T_A = L/vji, the Alfvenic crossing time, the solution is given by: 



1 — exp 



t 

r-R. 



\j ix,y,t)\ = \uph ix,y)\ (27rr) — 



1 — exp 

r-R. 



So that the magnetic energy Em and the ohmic dissipation J rate are given by 



1 

2,v 



L Uph 





1 — exp 1 















J 



TA 



1 — exp ( 



(4.81) 
(4.82) 



(4.83) 
(4.84) 



where Uph is the rms of Wph^ and with the rms of the boundary velocities u'^ and fixed 
to 1/2 we have Uph ~ 1. Equation (j4.84p implies that the heating rate for unitary volume 
grows quadratically in time for a time small compared to the resistive time ttj, while on 
this diffusive time scale the heating rate saturates to a value that is proportional to the 
value of the Reynolds number and the square of the axial Alfvenic velocity. 



4.7 Energy Balance 



As shown in § 12.1.21 the energy equation to which our system obeys is 

— = S-{J + n) (4.85) 



where 



E = \ l^d'x (ui + bi) , J = i l^d'xf, 17 = 1 l^d'xoj^ (4.86) 

S = (ldaS-n = +VA dau^-h±-VA dau°-b^, (4.87) 
/ Jz=L Jz=0 

E is the total energy, J and Vt are respectively the ohmic and viscous dissipation rates. S 
is the surface integral of the Poynting flux S (j2.79p 

S = B X (u X B) = B^u- (B • u)B (4.88) 

In reduced MHD (see § 12. ip the fields are decomposed in orthogonal and axial components 
B = t;_4e2 + b_|_, u = u^. The only dynamical fields are the orthogonal components of 
the magnetic and velocity fields b^ and u_|_, and in particular the velocity is supposed to 
have a vanishing axial component (at least in the first order contribution of the expansion 
described in ^ 12.11 in particular see (|2.32p ). So the only boundaries of the computational 
box that contribute to the surface integral S of the Poynting flux are the top {z = L, with 
velocity u^) and bottom {z = 0, with velocity u^) planes, while the contributions from the 
other planes cancel each other because of periodicity. This implies that the only term on 
the right hand side of equation (j4.88p that gives a contribution is the last one, because the 
other one has vanishing component along z. The axial component of S is then given by 

S, = S-e, = -VA (b^ • ux) (4.89) 
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and the total flux of energy from the boundaries S is given by (|4.87p . where the signs have 
been chosen so that S is positive when we have energy entering the system and negative 
when it leaves. 

One of the results of the anisotropic turbulence theory and an hypothesis of the 
reduced MHD expansion is that the spatial structures are elongated in the axial direction, 
because the cascade is strongly reduced in that direction. As shown in the linear analysis 
(see § 14. 2p the boundary conditions limit the velocity to its boundary value, while the mag- 
netic field grows to a higher value (see equations ( [^.5'7| j- |^.^^[ jj, so that in both the linear 
and nonlinear stages the magnetic field is the dominant field and its variation along the z 
direction is weak. 

We can so approximate the magnetic field b_L as uniform along the axial direction 
in the flux integral (|4.88p : 

S^VA J da (u^ - u°) • (4.90) 

and indicating with \iph = — u'^, with u'^ ^ but / da |u^|^ = Jda = 1/2, we 

have 

S r^v^J dauph ■ b_L (4.91) 

2D spatial periodicity in the orthogonal planes allows us to expand the velocity and magnetic 
fields in Fourier series, e.g. 

27r 

u (x, y) = Y, ^r,s e^'"'---^, where k^,, = — (r, s, 0) (4.92) 

r,s 

Then using this expansion in (|4.9ip we have 

Sr-VAYl (%^)r.,. • ff d:Edy bxe''^'--" = fvAY^ {uph),^^ ■ b_,,_, (4.93) 
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This equation shows that energy flux entering the system is proportional to the sum of 
the scalar product of the Fourier amplitudes of the velocity and magnetic fields at the 
same absolute wavenumber k,. ,,. From observations we know that most of the energy of 
photospheric motions is at the characteristic scale of a granule while shorter wavelengths 
have a reduced power level. Supposing that a turbulent cascade takes place, this implies that 
the Fourier components of the magnetic field hr^s are smaller at higher wavenumbers. In 
this way both the velocity and magnetic field amplitudes entering in the sum \4-93 ) decrease 



at high wavenumbers, so that the contributions to the flux S from scales smaller than the 
injection one are noticeably smaller and we neglect them. 

We now use the linear solution (j4.5ip - (|4.56p for the two-sided problem with velocity 
patterns u'' and u^, respectively in the boundary planes z = and z = L. Averaging over 
time-scales bigger than the Alfvenic crossing time, the magnetic field for this solution can 
be approximated as 

h± {x, y, z,t) ^— [u^ (x, y) - u° {x, y)] ~ — Uph (4.94) 

and from (j4.87p we can write for the injection energy rate 

S r^VA [da (u^ - u°) -h^^VA Ida |u^ - u°|^ — ~ ^^u^u^^ — (4.95) 

so that during the linear stage the injection energy rate is always positive. 
Taking into account diffusion, from (j4.80p - (j4.8ip we have that 



bj_ ~ Uph — 

TA 



t 

1 — exp 

' T-Jl 



so that from (jirsTl) 

t 



S (. VAU h 

TA 



1 — exp 



Tn 



(4.96) 



(4.97) 
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Figure 4.1: Poynting flux integral S, representing the energy rate entering the box, and the 
dissipation rate J. The parameters used are £ = 1, L = 10, u^h = ^-, va = 200 and TZ = 400 
to which correspond an Alfvenic crossing time t_4 = Ljvj^ = 5 • 10"^ and a diffusive time 
= 400/(47r • 4)2 ~ 6.3 • 10~\ and a saturation level Jsat = Ssat ~ 2.5 • 10^. 



For times smaller than the diffusive time t t-ji the injection of energy rate S is alway 
bigger than the dissipation rate J (|4.84p 

Sr^fvAul^—, J ^evAul^, — — (4.98) 

until for t ^ t-ji they both saturate to the value 

Jsat = Ssat = VAulj,^ (4.99) 

TA 

as shown in Figure 14. 1[ 
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Chapter 5 

Numerical Simulations 

Since Parker |67l [72] proposed his nano-flare scenario for coronal heating a num- 
ber of numerical experiments have been carried out to investigate this idea, with par- 
ticular emphasis on its possible relationship with the power-law distribution of observed 
emission events at optical, ultraviolet and x-ray wavelengths of the quiet solar corona. Mi- 
kic et al. [56], Longcope & Sudan [53j and Hendrix & Van Hoven [32] carried out simulations 
using a 3D "straightened out" loop and imposing photospheric shearing given by alternate 
direction flow patterns. They showed that a complex coronal magnetic field results from 
the photospheric field line tangling and although the field does not — strictly speaking — 
evolve through a sequence of static force-free equilibrium states, magnetic energy nonethe- 
less tends to dominate the kinetic energy. In this system, which may be thought of as 
evolving in a special regime of MHD turbulence, the field is restructured into current sheets 
elongated along the axial direction. 

The numerical investigation of the Parker scenario has been challenging since the 
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beginning. The reasons are the large aspect ratio of a coronal loop, and because the two 
different time-scales — along the axial direction (the Alfvenic crossing time ~ L/vji^) and 
the perpendicular direction (the non-linear time T^r/, ~ ^/uph), with t_4 ^ Tjvl — together 
require a large number of grid points. High-resolution simulations are also a necessity when 
investigating turbulence since this takes place at high Reynolds numbers. Turbulent systems 
are in fact characterized by a transition to the turbulent state which is not detected at small 
Reynolds numbers. A limitation of the first simulations carried out during the late 1980s 
and the early 1990s was their low resolution. While the results of these simulations mostly 
agreed among themselves, the interpretations that were given did not, and there was not a 
clear evidence in favor of one of them. 

To investigate the turbulent nature of this model with higher resolution, and to 
carry out long-time simulations necessary to define the statistics of heating events, Ein- 
audi et al. |22J first carried out 2D numerical simulations of incompressible magnetohydro- 
dynamics (MHD) turbulence using a random large scale magnetic forcing function to mimic 
the forcing exerted in three dimensions by the photosphere. The simulations displayed 
flux tube behaviour similar to the first 3D simulations, confirmed also in 2D simulations 
by Dmitruk et al. [16]. These 2D simulations were extended to longer times by Geor- 
goulis et al. [32], and showed for the first time how the magnetically dominated turbulence 
in the 2D system displays bursts in the dissipation (intermittency) which follow a power 
law behaviour in total energy, peak dissipation, and duration with indices similar to those 
determined observationally in X-rays. 

Recently Gudiksen &: Nordlund |39[ I40| have performed numerical simulations of 
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the Parker scenario, modeling a large part of the solar corona with a more realistic geometry. 
While this approach has advantages when investigating the coronal loop dynamics within 
its neighboring coronal region, modeling a larger part of the solar corona drastically reduces 
the number of points occupied by the coronal loops. Thus, these simulations have not been 
able to shed further light on the physical mechanism responsible for the coronal heating. 

In order to further investigate the underlying physical processes at work, we return 
to the simpler cartesian 3D model, performing numerical simulations with the highest reso- 
lution and longest duration to date. In section ^ 15.11 we describe the numerical simulations 
performed with a vortex-like velocity pattern and in § 15.21 those performed with a sheared 
velocity pattern. In § [6] we analyze in the framework of anisotropic turbulence theory the 
results of the simulations. 

5.1 Vortex-like Velocity Pattern 

In this paragraph we show the results of 3D numerical simulations modeling a 
coronal layer driven by a forcing velocity pattern constant in time. On the bottom and 
top planes (respectively z = and z = L) we impose two independent velocity forcings, 
resulting from the linear combination of large-scale eddies. The velocity potential at each 
boundary is given by 

fix, y) = ^ V - — ^==^ sin [2tt {kx + ly) + 27r ^ki] (5.1) 

where < x,^ < 1 and < z < L, with L = 10 (corresponding aspect ratio equal to 10). 

1 /2 

The wave number values {k, I) and (m, n) used are all those in the range 3 < (/c^ + /^) < 4, 
so that the injection wavenumber is km ~ 4 and the corresponding injection scale Im ~ 1/4. 
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n No. 
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VA 




n, h 


7^, 7^4 
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800 


oo 
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2 
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200 
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400 


oo 
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3 


V 


200 


128 X 128 X 50 


n 


200 


10 


2172 


4 


V 


50 


512 X 512 X 200 


h 


3 • 1020 


10 
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5 


V 


200 


512 X 512 X 200 
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10 


453 


6 


V 


400 


512 X 512 X 200 


h 


1020 


10 


77 


7 


V 
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512 X 512 X 200 


h 


1019 


10 


502 


8 
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200 


512 X 512 X 200 


n 


800 


oo 


551 
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200 


256 X 256 X 100 


n 


400 


oo 


4123 


10 


s 


200 


128 X 128 X 50 


n 


200 


oo 


2356 



Table 5.1: Summary of the simulation runs, v or s indicate the kind of boundary velocity 
forcing pattern: vortex or shear, va is the axial Alfvenic velocity value and rix x Uy x is 
number of points for the numerical grid, n, h indicate normal or hyperdiffusion. IZ or IZ4, 
indicates respectively the Reynolds number value and the hyperdiffusion coefficient, while 
7^11 is the Reynolds number along the axial direction where always normal diffusion is used. 
The duration of the simulation tmax/TA is given in Alfvenic crossing time unit = L/v_a. 



Oki and S^ki are two sets of random numbers whose values range between and 1 



0<aki,Cki<l, (5.2) 



and are independently chosen for the two boundary surfaces. The velocity field follows 
from equation ()4.4p . and the normalization in equation ()5.ip has been chosen so that the 
rms value of the corresponding velocity is l/\/2, i.e. 

dxdy {ul + uD = ^ (5.3) 

The initial condition is given by the strong uniform axial field specified by the value of the 
corresponding Alfvenic speed va (|4.3p . We have performed simulations with different values 
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(a) (b) 

Figure 5.1: Boundary conditions for our numerical simulations are specified imposing the 
velocity potentials ip^{x,y) in the bottom plane z = and ip^{x,y) in the top plane z = L. 
These result from the linear combination of large-scale eddies (wavelength ~ 1/4) with 
random amplitudes ()5.ip . normalized so that the velocity rms is I/V2. Fi gure shows the 
contours of the velocity potential in z = (a) and z = L (b) for the numerical simulation 
with v_/[ = 200 and TZ = 800. In lighter vortices the velocity field is directed anti-clockwise 
while in darker vortices it is directed clockwise. 

of VA- 

5.1.1 Simulations with v_a = 200 

We start showing the results of a numerical simulation performed with = 200, 
a numerical grid with UxXnyXriz = 512 x 512 x 200 points and a Reynolds number TZ = 800 
and 7^11 = 10. In Figure [5TT] we show the two boundary velocity patterns (equations (j5.ip ) 
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Figure 5.2: High-resolution simulation with = 200, 512x512x200 grid points and TZ = 
800, 7^11 = 10. (a) Magnetic (Em) and kinetic (Ek) energies as a function of time (r4 = 
L/v_A is the axial Alfvenic crossing time). Kinetic energy is always a small fraction of total 
energy, (b) Poynting flux S dynamically balances the Ohmic (J) and viscous (O) dissipation. 
Similar to kinetic energy dynamics, also enstrophy {Vt) is always a small fraction of total 
dissipation, confirming that the system is magnetically dominated. 

for this simulation which result from a specific random choice of the amplitudes ()5.2p . The 
total duration is roughly 500 axial Alfvenic crossing times {ta = L/vji,). To remark how 
challenging these kind of simulations are this simulation has used 52.428.800 grid points and 
814.215 time-steps which, using a 3^'"^ order runge-kutta, correspond to 2.442.645 substeps 
in which the derivatives have been computed in all the grid points. Fig. 15. 2al shows the total 
magnetic {Em) and kinetic {Ek) energies in the loop as a function of time, and Fig. I5.2bl 
shows both ohmic (J) and viscous ($7) dissipations and the Poynting flux S (the energy 
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that for unit time is injected into the system). Ah these quantities at first grow following 
roughly the linear behaviour described by the equations (I4.5ip -( j4.56p . until time t ~ 6x4 
when nonlinearity sets in. The system results magnetically dominated both for energy 
(magnetic energy E^j is ~ 35 times bigger than kinetic energy Ek) and dissipation (ohmic 
dissipation rate J is ~ 6.5 times viscous dissipation Q). 

In the fully nonlinear regime a statistically steady state is reached, in which the 
Poynting flux S, i.e. the energy which is continuously injected for unit time into the system 
at the boundaries by the field-line tangling due to the photospheric forcing, balances the 
energy which is dissipated for unity of time. On the other hand, as previously stated, the 
Poynting flux, depends not only on the forcing velocity imposed at the boundaries, but 
also on the value of the magnetic field generated inside the system. This is the reason for 
which this system must be considered a self-organized system. In fact the balance between 
Poynting flux and dissipation is reached by letting the magnetic field grow at the opportune 
value. 

In the absence of dynamical evolution, magnetic energy and ohmic dissipation 
would follow the curves (|4.83|) - (|4.84|) saturating on a time-scale given by t-ji = 7^/(27rn)^ ~ 
25r_4 respectively to the values Em ~ 3200 and J ~ 5000, well beyond the levels reached 
in the simulation. 

In order to verify the temporal stability of the steady state found in the fully 
nonlinear regime we have performed other numerical simulations, always with the value 
vj( = 200 for the axial Alfvenic speed, but with lower resolutions and Reynolds numbers 
to reach longer durations, namely x Uy x riz = 256 x 256 x 100 with TZ = 400 and 
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Figure 5.3: Two simulations with Alfvenic velocity = 200 but performed with lower 
resolutions to reach longer durations. In blue lines are shown physical quantities for a 
simulation with TZ = 200 and 128 x 128 x 50 grid points, while in green those for TZ = 400 
and 256 x 256 x 100. Both simulations do not have parallel diffusion T^y = oo. (a) Magnetic 
{Em) and kinetic (Ek) energies as a function of time. Kinetic energy is always a small 
fraction of total energy, (b) Ohmic (J) and viscous (fi) dissipation rates. Similar to kinetic 
energy dynamics, (O) is always a small fraction of total dissipation. Not shown in Figure 
the Poynting flux S always balances the dissipation rate. 



Ux X ny X riz = 128 x 128 x 50 with TZ = 200. In this earlier simulations no diffusion along 
the axial direction was used, i.e. T^y = oo. Time durations are respectively At ~ 1100 r_4 
and At ~ 2200 r^. In Fig. 15.31 we show the magnetic and kinetic energies, and the ohmic 
and viscous dissipation rates as a function of time for these simulations. The results of 
the previous simulation, which was carried out with a higher resolution, but for a shorter 
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duration, are fully confirmed, and in particular in the nonlinear regime a steady state is 
reached. In particular the integral quantities, i.e. total energies and dissipation rates, are 
intermittent in time. 

Fig. I5.3bl shows another very interesting result. Ohmic and viscous dissipation 
rates for the simulations with TZ = 200 and 400 roughly overlap. Having a dissipation 
rate independent of the Reynolds number is a property of turbulent systems where the 
injection, transfer and dissipation rates are equal (see eq. (|2.105p ) and independent of the 
Reynolds numbers beyond a threshold. The value of the threshold is determined by the 
diffusion time at the large scale, which must be larger than the nonlinear time-scale. The 
diffusion time at the scale A is r^) ~ T^A^, so that at the injection scale £j„ ~ 1/4 for our 
lowest resolution simulation with TZ = 200 the diffusion time is td ~ 250r_4 which we can 
suppose to be beyond the nonlinear timescale. On the other hand both this time-scales 
depend quadratically on A, so that with our lowest resolution we might be just beyond the 
threshold. 

We use hyperdiffusion (see (|3.3p ) to eliminate the diffusive effects at the large 
scales with a dissipativity exponent a = 4. In this case the diffusive time at the scale A 
is given by Tq ~ TZa A^". Numerically we require the diffusion time at the resolution scale 
Amin. = 1/-^) where N is the number of grid points, to be the same for both normal and 
hyperdiffusion, i.e. 



So for a grid with 512 points — which has required a Reynolds number TZ = 800 with 
normal diffusion — and a dissipativity a = 4, the required hyperdiffusion gives ~ 10^^. 



iV2 



a 




(5.4) 
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From equation (j2.73p we find that for a generic dissipativity a, including also a diffusive 
term along the axial direction z with Reynolds number TZu , dissipation rate is given by 



In general because of the anisotropy of the system Py <C 2?", for example in the high- 
resolution simulation shown in Figure 15.21 we have P|| ~ IQ-^VI even though the values 
of the Reynolds numbers used are T^y = 10 <C 7^ = 800. For normal diffusion a = 1 the 
dissipation rate (see equations (j2.73p - (|2.78p ) is approximated by P ~ ~ J + ri, with the 
ohmic and viscous dissipation rates as defined in (|2.84p - (|2.85p . 

We have then performed a numerical simulation using hyper diffusion, with dissi- 
pativity a = 4, Ux X Uy X Hz = 512 x 512 x 200 grid points. After a few tests we have 
used i?4 ~ 10^^ and a Reynolds number for the parallel diffusion T^y = 10. This simulation 
confirms the results with ordinary diffusion. In Fig. 15.41 we plot the dissipation rates for 
the 4 simulations with 3 different Reynolds numbers and hyperdiffusion. The simulations 
with TZ = 200 and 400 do not use parallel diffusion (T^y = oo) while those with TZ = 800 
and hyperdiffusion have 7?.y = 10. The dissipation rates roughly overlap, in particular those 
performed with the same numerical methods, i.e. those which implement or not parallel 
diffusion. We also note that these four simulations do not differ solely for the value of the 
Reynolds numbers, the number of grid points, and the implementation of parallel diffusion: 
in each simulation the coefficients that define the spatial shape of the forcing boundary 
velocities (see eqs. (j5.ip - (j5.2p ) are chosen randomly, and are always different. Thus, the 
four simulations shown in Fig. I5.4al refer to four different velocity patterns, although all 




(5.5) 
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(a) (b) 

Figure 5.4: (a) Dissipation rates as a function of time for four numerical simulations with 
different Reynolds numbers and hyperdiffusion with dissipativity a = 4. Their close values 
over a very large range of diffusion times at the large scales show that dissipation is inde- 
pendent of the Reynolds numbers, a characteristic of turbulent systems, (b) Close-up of 
dissipation rates. With increasing Reynolds number shorter temporal variations develop. 

of them have an rms value of (see (15. 3p ). The small differences are probably due 

to this difference more than to the different Reynolds numbers. To further investigate the 
independence of the dissipation rate on Reynolds numbers it will be necessary to perform 
these numerical simulations using the same spatial pattern for the boundary velocities and 
parallel diffusion for all of them. The diffusion time at the large scales has a tremendous 
increase between TZ = 200 and TZ^ = 10^^, going from td ~ 250 r_4 up to T£)_4 ~ 3 • 10^^ r_4, 
so that the results shown in Figure 15.41 demonstrate that there is at most a very small 
variation in the dissipation rates once the diffusion time at the large scales is beyond a 
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Figure 5.5: Numerical energy conservation error jA/S"! (see eq. ()5.6p ) as a function of time, 
threshold. 

Figure I5.4bl shows a close-up of the dissipation rates as a function of time at 
different Reynolds numbers and with hyperdiffusion. At increasing higher Reynolds numbers 
the function clearly exhibits more fine structures at smaller time-scales. This is a typical 
feature of a transition to turbulence. 

In order to check the validity of these results it is important that the energy equa- 
tion ()4.85p is numerically well satisfied. Although this equation is not directly integrated 
by the numerical code, it follows from the equations that are integrated, so that there could 
be a numerical discrepancy. Indicating with A the quantity 

r) F' 

A = — -[S-V], (5.6) 

which analytically should be equal to 0, we define as the numerical energy conservation 
error the absolute value of the normalized quantity A/5. In Figure [S3] we plot |A/5| as 
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a function of time for the numerical simulations just described. The value of |A/5| is at 
most 2.5%, with an average of 1.8% for n = 200 decreasing to 0.15% for 7^4 = 10^^ while 
for TZ = 400 and 800 the average errors are respectively 0.9% and 0.4%. Hence the energy 
conservation equation (|4.85p is numerically very well satisfied. 

We now briefly examine the 3D structure of the physical fields. One of the hy- 
pothesis of reduced MHD is that during time evolution the orthogonal magnetic field is 
small compared to the axial component. The orthogonal field component always fluctuates 
around a value which is roughly the 3% of the the axial field (see Figure I5.2ap thus verifying 
this hypothesis. In Fig. 15.6115.71 we show the magnetic field lines, a view from the side and 
from the top of the computational box respectively. For an improved visualization the box 
has been rescaled, but it should be noted that the computational box is ten times longer in 
the axial direction {z) than in the orthogonal plane. Consequently it can be verified that 
the magnetic field lines are only slightly bent. In both figures we show the field lines of the 
orthogonal magnetic field in the mid-plane. The structure of the orthogonal magnetic field 
is almost invariant in the axial direction and it is structured in magnetic islands. In partic- 
ular, no boundary layer is present. The lines of the total magnetic field which bend most 
are those which happen to be at the outskirts of the magnetic islands, where the orthogonal 
magnetic field is enhanced. Current sheets develop, extended along the axial direction, as 
shown in Fig. 15. 8115. 91 and in these thin sheets the energy flux is finally dissipated. 

5.1.2 Spectral Properties 

In this paragraph we discuss the dynamics of our system using its spectral quan- 
tities in more detail. We show how the current sheets shown in Figures I5.8ti5.9l are formed 
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Figure 5.6: Field lines of the total magnetic field (orthogonal plus axial) at time t = 18.5 r_4 
for the high-resolution numerical simulation with Vj{ = 200, 512 x 512 x 200 grid points, 
TZ = 800 and T^y = 10. Mid-plane: field lines of the orthogonal magnetic field. The 
orthogonal magnetic field magnitude fluctuates around a value which is roughly the 3% of 
the axial component, well within the reduced MHD ordering. This is reflected in the slight 
bending of the magnetic field lines. For an improved visualization the box size has been 
rescaled. But the axial length of the box is ten times longer than the orthogonal one. The 
resize of the box artificially enhances the field line bending. The orthogonal magnetic fields 
is structured in magnetic islands, and is mostly homogeneous in the axial direction (z). 
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Figure 5.7: A top view of the magnetic field lines shown in Figure 15.61 Mid-plane: field 
lines of the orthogonal magnetic field. The orthogonal magnetic field is structured in mag- 
netic islands, and is mostly homogeneous in the axial direction (z). Because of the quasi- 
homogeneity of the orthogonal field, the field lines of the total magnetic field are more bent 
when they are at the outskirts of magnetic islands, while at different heights they always 
have the orthogonal component of the magnetic field roughly in the same direction. 
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Figure 5.8: Isosurfaces of the squared current j^: side view. Two isosurfaces of the squared 
current at time t ~ 18.5 for the numerical simulation with = 200, 512x512x200 grid 
points and Reynolds numbers TZ = 800, T^y = 10 is represented. The isosurface at the value 
j'^ = 2.8- 10^ is represented in partially transparent yellow, while red displays the isosurface 
with = 8 • 10^, well below the value of the maximum of the squared current that at this 
time is j'^ = 8.4 • 10^. N.B.: The red isosurface is always nested inside the yellow one, and 
appears pink in the figure. 
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Figure 5.9: Top view of the isosurfaces of the squared current shown in Figure 15.81 using 
the same color display. The isosurfaces are extended along the axial direction, and the 
corresponding filling factor is small. 
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and study the long-time dynamics behaviour. 

In turbulence the fundamental physical fields are the Elsasser variables = 
u_L lb b_L. The analysis of the weak and strong anisotropic turbulence presented in § 12.2.41 
has been carried out assuming weak velocity-magnetic-field correlation at all scales 6z^ ~ 
6z^ ~ 5zx. This property basically follows from the symmetry in of both the reduced 
MHD equations p.53p - ()2.56p and the boundary conditions used p.60p - ()2.6ip . as shown in 
§ 16.11 and will be used in all § E] to perform our nonlinear analysis. The Elsasser variables 
are associated with the corresponding energies 

E^ = - [ d^x (z±)' (5.7) 
2 Jv 

related to the kinetic and magnetic energies Ek, Em and to the cross helicity H'-^ 

H^' = [ d\ • bx, (5.8) 
Jv 

by 

= Ek + Em ± (5.9) 

so that the low-correlation of the fields requires < Ek + Em- Fi gure I5.10al shows this 
quantity as a function of time for the simulation with vji^ = 200, TZ = 800, TZ\\ = 10. The 
cross helicity fluctuates around zero, with a variation of at most 4% with respect to total 
energy and an average of 0.54%. Alternatively for a magnetically or kinetically dominated 
systems, in which one of the two energies exceeds the other, the correlation function is 
defined as 

P = m = — ^-10 

{^EkEm)^'^ (/d3x • /d3x b^ 
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(a) (b) 

Figure 5.10: (a) Ratio between cross helicity and total energy ii^ as a function of time, 
(b) Velocity-magnetic field correlation p'^ as a function of time, showing that the system is 
weakly correlated. Both quantities are calculated for run 1. 

This quantity is shown as a function of time in Figure IB. lObI It fluctuates around zero with 
an average of 0.017 and a maximum of ~ 0.1. Hence the velocity-magnetic-field correlation 
is very small when the whole computational box is considered. We present a more detailed 
correlation analysis later in this paragraph. In Figures [5.111 we show the spectra for run 1. 
Unless otherwise specified the spectra are one-dimensional, averaged in the axial direction. 
We define the Fourier transform of a function f{x,y) as 

e 

/(^) = ^/(fc)e'^-, f{k) = ^ll dxdy/(a;)e-^-, (5.11) 
k 

where 

27r 



Ill 



The discrete form of Parseval's theorem gives: 

dxdy |/(^)p = £2^|/(fc) ' (5.12) 



k 





For example considering the energy spectra 

£ 

E^ = ^ J'^dz jjdxdy{z^f = \ y"''dz^2^|z±|'(k,z) (5.13) 

° ° k 

Even if k is discrete in our case, it is useful to use a continuum formalism through the 
substitution 

(t) 5 ^ / ^^-^^^ 

In this way we can write the energies as: 

E^ = ^£dz (^^^ JJdk^dky |z±|'(k,z) = ^1^^ £dz2TTjdkk |z±|'(fc,z) 

(5.15) 

and finally integrating along z we have E^ = Jdk E^ , and in the discrete case 

E^ = Y,Et, n = l,2,... (5.16) 

n 

In Figure [S.llal we show the spectrum E^ for run 1. In different colors the spectrum at 
different times is displayed, while the black line shows the time average of the spectrum 
performed on a total duration of ~ 500 Alfvenic crossing times. A well-developed spectrum 
is formed with a peak at n = 4, which is the injection scale. This is a clear proof that this 
system is turbulent and that an energy cascade exists. Figure [5. llbl shows the time averages 
spectra of the kinetic, magnetic and E^ energies. As expected E^ almost completely 
overlap, while they are barely distinguishable from the magnetic energy spectrum. The 
kinetic energy spectrum has noticeably smaller values, and in particular the mode n = 4 
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Figure 5.11: (a) E~ energy spectra at different times are shown in color. Their time average, 
performed over a duration of ~ 500 Alfvenic crossing times, is shown in black, (b) Time 
averaged spectra for magnetic (Em), kinetic [Ek) and E^ energies. 



at the injection scale is fixed by the boundary forcing velocity as we show later in this 
paragraph. 

Figure [5.121 shows a close-up of the time-averaged E~ spectrum. Even if run 1, 
for which this spectrum has been calculated, used a grid with 512 x 512 x 200 points and 
a Reynolds number TZ = 800 an inertial range is barely formed and diffusion affects the 
spectrum at least up to n ~ 10. On the other hand energy is injected at the scale n ~ 4 at 
small wavenumbers, so in this kind of simulation it is very difficult to obtain a well-developed 
inertial range because very few points are left. Two lines corresponding to energy spectra 
kj_ and are drawn for reference. In order to study the properties of the inertial 
range we have used hyper diffusion, that we will describe in a moment. But now we want to 



113 



ioa.oooF ' ' " • ' ■ • I 




10 1D0 

n {wavenumber) 



Figure 5.12: Close-up of time-averaged E energy spectrum. An inertial range is barely 
formed, even using 512 x 512 x 200 grid points corresponding to a Reynolds number TZ = 800. 

analyze how the current sheets elongated along the axial direction shown in Figures [5.8115.91 
are generated. 

In Figures I5.13al and I5.13cl we show the first 7 modes and mode n = 1 1 for the 
magnetic and kinetic energy spectra, for the first 20 Alfvenic crossing times t_4. In Fig- 
ures 15.141 and 15.151 we show at selected times, covering the same time-interval, the axial 
component of the current j and vorticity oj in the midplane z = 5. The modes at the injec- 
tion scale n = 4 follow the linear dynamics described by equations (|4.5ip - (|4.56p until time 
~ 6t^. In particular the magnetic field and the current at time t ~ 0.63r_4 — i.e. just after 
time t = 0.5 r_4 when the two counterpropagating front- waves reach the midplane — are 
mapping the photospheric velocities — u'' as predicted by our linear analys. Figures [5. lal 
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Figure 5.13: Magnetic (a)-(b) and kinetic (c)-(d) energy modes as a function of time, (a)- 
(c) are a close-up, showing the first 20 Alfvenic crossing times r_4, of the respective figures 
(b)-(d) which show the dynamic of this modes for the whole duration of run 1 ~ 500 r4. 
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Figure 5.14: Axial component of the current j (in color) and field-lines of the orthogonal 
magnetic field in the mid-plane (z = 0), at selected times covering the time interval shown 
in Figures I5.13al and I5.13cl 
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(c) (d) 

Figure 5.15: Axial component of the vorticity cj (in color) and field-lines of the orthogonal 
magnetic field in the mid-plane {z = 0), at selected times covering the time interval shown 



in Figures I5.13al and I5.13cl 
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(a) (b) 

Figure 5.16: (a) Contour lines of the linear combination of the boundary velocity fields 
_ -qO ^Ij-j Qui-i-gnt density j squared averaged in the orthogonal x — y planes as a 
function of z at selected times. Same averages for the orthogonal magnetic and velocity 
fields. 



and lS.lbl show the two photospheric velocity patterns separately, while for comparison with 
Figure I5.14al we show the contour of the linear combination — u*^ in Figure I5.16al 

In particular the magnetic energy mode n = 4 grows roughly quadratically in time, 
while the n = 4 kinetic energy mode exhibits a sawtooth structure with a period r4. The 
dynamics are described by a cascade towards smaller scales. No instability of any kind is 
detected. In physical space this cascade corresponds to the formation of small scales which 
are organized in vortex-current sheets. This picture remains unaltered throughout the rest 
of the nonlinear stage. Energy is injected at the large scales (n ~ 4) by photospheric motions 
and an energy cascade develops, which transports this energy from the large to the small 
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scales where it is finally dissipated. Figures [5.14cl and I5.15cl show that during the nonlinear 
stage (t ~ 18.5 r_4) the magnetic field exhibits a complex topology and that distorted 
current sheets are almost always present together with quadrupolar vorticity structures. 
These current sheets are subsequently subject to tearing, which leads to suppose that the 
dissipation mechanism is nonlinear magnetic reconnection in a turbulent environment. 

An important question is what is the dissipation rate and how much it is influenced 
by the magnetic reconnection rate. Longcope &: Sudan [53] suppose that the reconnection 
proceeds at a rate independent from the turbulent state of the system, which they completely 
neglect. 

A similar behavior, i.e. microcurrent sheets formed by a turbulent cascade and 
then subject to tearing, has been detected in decaying 2D turbulence e.g. by Politano 
et al. |7B] . Biskamp &; Welter and Biskamp &: Schwarz |S| have also studied its influence 
on dissipation rates. In particular Biskamp & Schwarz [S] have performed the highest 
resolution 2D MHD simulations finding that the tearing mode occasionally occurs but its 
influence is rather weak. This is because this reconnection is taking place in a turbulent 
environment, so that the current sheets dynamics is influenced by the larger-scale dynamics. 

Our numerical simulations leads us to suppose that whether or not the reconnection 
rate is influenced by the turbulent dynamics the tearings happen continuously, so that the 
energy flux trasported by the cascade is dissipated in many "microbursts" , so that the 
turbulent energy flux is dissipated in a statistically steady fashion. 

In particular, the fact that the spectrum which is developed (see Figure [5. 11 1) does 
not show any anomalous structure at smaller scales means that the rate at which energy is 
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dissipated by magnetic reconnection is that at which energy flows along the inertial range. 
This scenario is reinforced by the fact that when we use hyperdiffusion the spectra that are 
found are consistent with those expected theoretically for weak anisotropic turbulence, as 
shown in the following paragraphs. 

The only difference detected at longer times is the growth of the n = 1 and 
n = 2 modes for the magnetic energy (see Figures I5.13bp , which show a tendency towards 
an inverse cascade, even though a complete inverse cascade never fully develops. This is a 
difference with the corresponding 2D simulations (Einaudi & Velli [21]), where at long times 
a complete inverse cascade leads to the coalescence of the magnetic islands. The kinetic 
energy mode at the injection scale n = 4 fluctuates around an equilibrium value which is 
the average of its value during the linear regime, and indicates that the dynamics for the 
modes at n = 4 do not depart dramatically from the linear behavior ()4.5ip - (|4.56p . 

The inhibition of an inverse cascade is probably due to the stiffness of the magnetic 
field lines, which are only slightly bent along the axial direction. In fact the magnetic field 
line tension associated with these field lines inhibits the coalescence of the magnetic islands 
that form in the orthogonal planes. In Figure [5. 171 we plot the first 4 modes for the magnetic 
energy with two different values of the axial Alfvenic velocity, respectively = 50 and 
vj_ = 1000. At VjX = 50 the field line tension is smaller than for vj( = 1000, and in fact in 
the first case an inverse cascade is clearly present (Figure I5.17ap while in the second it is 
not (Figure [5. 17bp . 

We show in Figure I5.16bl some physical quantities averaged in the x — y planes 
as a function of the axial direction z. Both the magnetic and velocity fields are elongated 
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(a) VA = 50 (b) VA = 1000 



Figure 5.17: Magnetic energy modes n = 1, 2, 3, 4 (n = 4 is the injection wavenumber) as 
a function of time for run 4 (a) and run 7 (b), respectively with v_a = 50 and vji, = 1000, 
showing that a bigger axial magnetic field inhibits through magnetic field line tension an 
inverse cascade. 



along the axial direction, as is the current. 

Finally in Figure I5.18al we show the magnetic helicity density integrated in the 
X — y planes as a function of z at selected times 

h'^ {z,t) = J J dxdyu^-h± (5.17) 



The figure shows an almost linear behavior of /i*^ as a function of z, confirming the homo- 
geneity of this system in the axial direction, /i*^ is the axial component of the Poynting 

vector S divided by the Alfvenic velocity v_a (see eq. (|4.89p ). i.e. 

e 

h'^' {z,t) = -^ jj dxdyS-e^ (5.18) 
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(a) (b) 

Figure 5.18: (a) Cross helicity density as a function of z at selected times, (b) Velocity- 
magnetic field correlation in the orthogonal planes as a function of z at selected times. 



so that energy is flowing from the boundaries into the computational box at a uniform 
rate. The energy equation (|4.85|) still holds for any slice of the computational box included 
between two planes z = const so for any two slices of equal volume the rate of injection of 
energy and the dissipation rate are the same. 

The velocity-magnetic field correlations in the x — y planes p'^y, which is closely 
related to the cross helicity density h'~^ , is defined as 

n , s f fn dx d-U U I • b I 
P^y {Z, t) = —y-^ (5.19) 

(//o dx dy ul ■ JJ^dx dy 
In Figure I5.18bl we show this quantity as a function of z at selected times. As expected 

the velocity and magnetic field are more correlated at the boundaries, where the velocity is 

imposed as a boundary condition, than in the interior of the computational box. 
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5.1.3 Weak Turbulence Spectra 

In order to resolve the inertial range we perfomed a series of high-resolution sim- 
ulations using hyper diffusion. The numerical grid was 512 x 512 x 200, and the axial 
Alfvenic velocities are vj\^ = 50, 200, 400, 1000. In Figure [5.191 we show the time-averaged 
E~ energy spectrum for these simulations. The time averages have been performed over 
~ 200, 450, 77 and 500 r^^ respectively. Comparing Figure [5.19bl with Figure [S.llal we 
notice that a well-developed inertial range is now present, in contrast to the analogous case 
with normal diffusion. In the framework of the weak anisotropic turbulence theory that we 
have briefly surveyed in § 12.2.41 the spectral slopes are expected to increase from — 5/3 for 
strong turbulence, i.e. for turbulent regimes in which the perturbation and the mean field are 
similar, to steeper exponents as long as we proceed towards weaker regimes of turbulence, 
i.e. turbulent regimes in which the ratio of the perturbation over the mean field decreases. 
For weak turbulence the exponents are a„ = —(3n — 5)/(n— 1), i.e. —2, —7/3, —5/2,.... To 
distinguish between spectral slopes is usually a difficult task because the differences among 
them are relatively small. In the case of weak turbulence at increasing n the difference 
between neighboring slopes an becomes very small. 

Consider Figure [5.19b|. which is the run with the same value of the Alfvenic velocity 
that we have extensively considered in the previous paragraphs vj^ = 200. This is the easiest 
case: while the k~'^ line fits well the inertial range, the two neighboring slopes k^^^^ and 
^-7/3 cigarly do not. 

A very well known problem with hyper diffusion is the so-called bottleneck effect 
(Falkovich |23] . Lithwick & Goldreich the spectrum on scales slightly larger than the 
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Figure 5.19: energy spectra for 4 simulations with = 50, 200, 400, 1000. Hyper- 
diffusion with dissipativity a = 4 has been used for all the simulations with respectively 
7^4 = 3 • 10^'^, 1 • lO^'^, 1 • 10^^, 1 • 10^^. All the simulations use a numerical grid with 
512 X 512 X 200 grid points. 
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dissipation scales becomes flatter (in bilogarithmic scale a small bump appears). The energy, 
in effect, is backed up and the origin of this reflection is the rapid increase of diffusion at 
the small scales. Hyperdiffusion virtually removes diffusion from the large scales, but while 
normal diffusion slowly decreases the diffusion time towards the small scales, hyperdiffusion 
has a very steep decrease which increases with increasing dissipativity a. We have used a 
mild exponent a = 4 but the bottleneck effect is present visible, especially in the simulation 
with f_4 = 400 (Figure l5.19cp and 1000 (Figure l5.19dp . In Figure [5.19cl the spectrum is 
steeper than for vj^ = 200 and the slope that fits better the inertial range is while -5/3 
and —3 certainly do not. Consider now Figure [5.19dl where the spectrum for a simulation 
with the higher value v_a = 1000 is shown. Between —5/2 and —3, the first slope fits better 
in the inertial range but —7/3 and —2 clearly do not fit the inertial range. The case with 
the lowest value of the Alfvenic velocity v_a = 50 shown in Figure I5.19al The spectrum 
shown is an average over 196 Alfvenic crossing times. The —7/3 slope fails in the inertial 
range, but —2 and —5/3 seem to fit the first part and the last part of the inertial range 
respectively. Actually we could consider if in this case it does develop an inertial range. 
The slope change which develops around n ~ 10 is in a very advanced position towards 
the small wavenumbers respect to the others simulations to be the bottleneck effect. But 
we have no objective method to say if this is really a bottleneck effect. There is another 
possibility. As shown in § 12.2.41 all the weak turbulence scalings share the property that 
at smaller scales the turbulence gets stronger, and Goldreich &; Sridhar [34j have proposed 
the possibility of multiple inertial ranges, where starting, for example, from a —2 slope the 
system transitions at smaller scales to a stronger —5/3 regime. If it is already difficult to 
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Figure 5.20: The rms of the fields, which is roughly the same for the two fields and that 
we generically indicate with z, is shown as a function of time for the four simulations with 
vj\^ = 50, 200, 400 and 1000. In (a) the relative value z/vj\^ is shown as a function of time, 
while in (b) it is shown z. 

distinguish numerically one inertial range, it is certainly far more difficult to distinguish 
among two. Nevertheless in Figure I5.19al either we have a bottleneck effect or it is this 
transition, and there is no objective method to decide. It would be certainly interesting to 
perform a numerical simulation with the same parameters but doubling the resolution. 

As already shown, for this system ~ so that E'^ ~ E~ . Hence the rms of 
the fields are also roughly the same, and we now consider the integral quantity 



2 {Ek + Em) 1 /■ ,3 I ±,2 



(5.20) 



Figure I5.20al shows the relative ratio z/vji^ as a function of time while Figure I5.20bl show z 
for the four simulations considered in Figure 15.191 While the value of z increases at higher 
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Figure 5.21: Dissipation rates for the four simulation with vj^ = 50, 200, 400, 1000. 

values of w^^, the relative ratio z/v_a decreases. Hence the results shown in Figures [5 . 19115 . 2 J] 
show that the our system is described by a weak anisotropic turbulence regime which gets 
weaker for higher values of the axial Alfvenic velocity v_a. 

The physical mechanism by which the system dynamically choses a value z for 
different values of VJ!^ is self-organization. Both this mechanism and its consequences for 
coronal heating scalings are examined in § [6l We conclude this paragraph showing the 
dissipation rate as a function of time for the four simulations in Figure [5.21 1 the dissipation 
rate is higher for higher values of vji^. 

5.1.4 Strong Turbulence and Dissipation 

In previous paragraphs we have shown the spectra of many physical quantities 
always considering the orthogonal k±_ variable, and summing over the axial k^ direction. 
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(a) (b) 

Figure 5.22: (a) E~{nz) spectrum as a function of + 1. The small scale increase is due 
to aliasing, (b) E~ {n±,nz) as a function of + 1 at selected values of n±. Spectra are 
calculated for run 7 {R^ = 10^^^;^ = 1000) at time t ~ 500 r^. 

In our simulations we impose non-periodic boundary conditions, along the axial direction 
z, but nevertheless the system results almost periodic, i.e. there are not strong variations 
along the axial direction. It is then possible to perform a fourier transform of the physical 
quantities also along the axial direction. As the functions are not strictly periodic we have 
to expect some "aliasing" at the small scales. 

Consider the Fourier expansion (I5J5]1 for the energies. If instead of summing 
along z we perform the fourier transform we obtain the 2D analogue of the ID spectra (|5.16p . 
and we can now write 

^± = 5^ ^±(nx,n,), n^ = l,2,... n, = 0,l,2,... (5.21) 

where n± is the orthogonal wavenumber, and is the axial wavenumber for which also the 
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Figure 5.23: Logarithmic contour of the E~ {n±,nz) spectrum at time t ~ 500 r^^ for run7. 
At sufficiently high orthogonal wavenumbers n± > 20 the formation of smah scales along the 
axial direction becomes significant. The over-plotted line represents the curve oc nij^^ . 



Uz = component does not vanish. From the 2D spectra {n±,nz) we easily obtain the 
ID orthogonal and axial spectra E^{n±) and E^(nz) summing over the remaining variable, 
i.e. 

E^ (n^) = J2^^ (^^' ' (^-) = J2^^ ' (^-22) 

Uz n± 

In Figure [5.221 the spectra E~[nz) and E~ {n±,nz) at selected values of n± are 
shown. Figure I5.22al shows that already the first few modes have a sharp decrease in their 
values. Figure [5.22bl shows that at least up to n± ~ 20 the presence of axial (z) modes is 
negligible. 

In Figure !??] are shown the contours of the natural logarithm of the 2D E~ (n_|_, n^) 
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spectra. It is clearly shown that at small scales strong turbulence sets in, with the formation 
of small scales along the axial direction. 

5.2 Sheared Velocity Pattern 

In this section we consider the 3D evolution of a coronal layer driven by a time 
independent forcing pattern. A vanishing photospheric flow is imposed on the bottom plane 

{z = 0) 

u±{x,y,z = 0,t) =0, (5.23) 

while on the top plane {z = L) the velocity pattern is a stationary photospheric shear flow 
aligned along the y direction 

u±{x, y, z = L,t) = sin(27rnx + 1) (5-24) 

where n = 4, 0<x,y<l and < z < L, with an aspect ratio of 10, i.e. L = 10. 
The initial condition is given by the strong uniform axial field with very small amplitude 
noise in the perpendicular velocity and magnetic fields. This case is interesting because the 
growing magnetic field induced in the corona (see equations (j4.8ip - (j4.84p ) by such motions 
is an equilibrium field at every instant. In other words, any dynamical evolution beyond 
the slow increase of magnetic energy with time must be due, at first, to instabilities arising 
in the system. A similar configuration was used by [44J, who developed a semi-analytical 
quasi self-consistent "turbulent" model of coronal heating. The most interesting result from 
these simulations is that after the first energetic burst (due to a tearing instability), which 
releases the energy accumulated during the linear stage, the dynamics are turbulent once 
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the nonlinear stage sets in and are similar for all to those which develop with a vortex-like 
forcing velocity pattern. 

If there were no perturbation at all, at longer times diffusive terms would not be 
negligible even at large scales and the solution would be modified over diffusive time-scales. 
While diffusion will only slightly change the shape in z of the velocity field, it has a stronger 
effect on the magnetic field which would otherwise grow linearly in time. We describe briefly 
such effect. From which the following expressions for the total magnetic energy and ohmic 
dissipation are obtained: 



Em = - / d^a; 6^ = , 



(5.25) 



(5.26) 



(27rn)' 

At first, for a time short compared to the resistive time tti = TZ/ (27rn)^, both quantities 
grow quadratically in time while they saturate on the diffusive time-scale. If non-linearity 
did not set in, the heating rate would saturate at a value that is proportional to the Reynolds 
number and the square of the axial Alfvenic velocity. 

The results of a numerical simulation performed with UxXnyXUz = 512 x 512 x 200 
grid points, a Reynolds number TZ = 800 and carried on for roughly 550 axial Alfvenic 
crossing times {ta = L/va) ^-^e shown in Figs. 15.241 and 15.251 Fig. I5.24al shows the total 
magnetic and kinetic energies in the loop as a function of time. Fig. I5.24bl shows ohmic 
and viscous dissipation, and Fig. 15.251 shows the total dissipation and the Poynting flux. 
Magnetic energy grows quadratically at first, reaching a maximum at t ~ 30 while 
the kinetic energy is much smaller and remains limited at longer times. In the absence of 
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(a) (b) 

Figure 5.24: High-resolution simulation with 512x512x200 grid points and TZ = 800. Left: 
Magnetic (Em) and kinetic (Ek) energies as a function of time. In dotted line curve (j5.25p . 
which represents the linear saturation of the magnetic energy if there was solely the linear 
dynamics. At time t ~ SOr^ a reconnection event abruptly decreases the value of magnetic 
energy. Kinetic energy is always a small fraction of total energy. Until time t ~ 30 ta 
magnetic energy follows the linear saturation curve (j5.25p . Right: Ohmic (J) and viscous 
(17) dissipation as a function of time. Curve (|5.26p represents the linear saturation of the 
ohmic dissipation and is drawn in dotted line. At time t ^ 30 ta the reconnection event 
abruptly releases the energy previously stored in the orthogonal magnetic field. Similar 
to kinetic energy dynamics, also enstrophy (17) is always a small fraction of total dissipa- 
tion, confirming that the system is magnetically dominated. Until time t ~ 30 ta ohmic 
dissipation follows the linear saturation curve (|5.26p . 
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Figure 5.25: Total dissipation J + Q (ohmic plus viscous dissipations) and incoming Poynt- 
ing flux 5 as a function of time for a high-resolution simulation with a grid of 512x512x200 
points and R = 800. Sign convention is chosen so that flux is positive when energy goes into 
the system and negative when it goes out. After the reconnection event at t ^ 30 ta a sta- 
tistically steady state is reached where energy which is continuously injected for unity time 
into the system at the boundaries through the field-line tangling, due to the photospheric 
forcing, balances the energy which is dissipated. 
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dynamical evolution, magnetic energy and ohmic dissipation would follow the curves (|5.25p - 
()5.26p (as they actually do until t ~ 30 when nonlinearity sets in), saturating on the 
diffusive time-scale t-ji, and shown by the dashed lines in Fig. 15.241 What occurs at t ~ 30 
is that reconnection begins in the sheared induced field. Given the periodicity of the system, 
the actual instability is a multiple tearing mode, which grows faster than the classical tearing 
mode. Until t ~ 30 Tyi energy is stored in the large-scale magnetic field, while the kinetic 
energy is only a small fraction of the total energy. This energy is abruptly released in the 
reconnection event that gives rise to the first big burst in the heating rate that consequently 
lowers the magnetic energy. After this event a statistically steady state is reached. In this 
state the Poynting Flux balances the rate of energy dissipation, as shown in Fig. 15.251 On 
the other hand, as previously stated, the rate of energy injection into the system, i.e. the 
Poynting flux, depends on both the forcing velocity imposed at the boundaries and the 
magnetic field generated inside the system. This is why this system is self-organizing. The 
balance between Poynting flux input and dissipation rate is reached by letting the magnetic 
field grow at the "opportune" value, as it is shown in § [6j 

Now we examine briefly the 3D structure of the physical fields. One hypothesis of 
reduced MHD is that during time evolution the orthogonal magnetic field is always small 
compared to the axial component. After the first big burst, where the ratio of the rms of 
the orthogonal field over the axial field reaches the 6%, this component always fluctuates 
around a value which is roughly the 3%. Thus the reduced MHD requirement is satisfied. In 
Fig. 15.26115.271 we show the magnetic field lines, respectively a view from the side and from 
the top of the computational box. For an improved visualization the box has been rescaled. 
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Figure 5.26: Field lines of the total magnetic field, orthogonal plus axial. Mid-plane: field 
lines of the orthogonal magnetic field. The orthogonal magnetic field magnitude fluctuates 
around a value which is roughly the 3% of the axial component, well within the reduced 
MHD ordering. This is reflected in the slight bending of the magnetic field lines. For an 
improved visualization the box size has been rescaled. But the axial length of the box is 
ten times longer than the orthogonal one. The resize of the box artificially enhances the 
field line bending. The orthogonal magnetic fields is structured in magnetic islands, and is 
mostly homogeneous in the axial direction (z). 
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Figure 5.27: A view from the top of the field lines of the total magnetic field (orthogonal plus 
axial). Mid-plane: field lines of the orthogonal magnetic field. The orthogonal magnetic 
field is structured in magnetic islands, and is mostly homogeneous in the axial direction 
(z). Because of the quasi homogeneity of the orthogonal field, the field lines of the total 
magnetic field are more bent when they are at the outskirts of the magnetic islands, while at 
different heights they always have the orthogonal component of the magnetic field roughly 
in the same direction. 
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But it should be noted that the axial length is ten times the orthogonal length of the box. 
The magnetic field lines are only slightly bent. In both figures we show the field lines of 
the orthogonal magnetic field in the mid-plane. The structure of the orthogonal magnetic 
field is almost invariant in the axial direction and it is structured in magnetic islands. In 
particular no boundary layer is present. The field lines of the total magnetic field which 
bend more are those which happen to be inside the magnetic islands where the orthogonal 
magnetic field is enhanced. The magnetic island structure for the magnetic field gives rise 
to current sheets elongated along the axial direction, as shown in Fig. 15. 28115.291 where the 
energy fiux is finally dissipated. 

We have also performed other numerical simulations with different Reynolds num- 
bers and spatial resolutions. In Figs. 15.301 we show the results of a very long simulation 
with Ux >i ny X riz = 256 x 256 x 100 grid points, a Reynolds number TZ = 400 and carried 
on for more than 4000 r_4. The results of the previous simulation, which was carried out 
with a higher resolution, but for a shorter duration, are fully confirmed. In particular the 
longer interval allows us to confirm that after the first reconnection event a steady state is 
finally established. 

We have another very interesting result. In Fig. 15.311 we plot ohmic dissipation 
for 3 numerical simulations carried out with 3 different Reynolds numbers, TZ = 200, 400 
and 800. As the Figure clearly shows they overlap. This supports our hypothesis that the 
dynamics are turbulent and that the dissipation is independent of the Reynolds number. 
In Fig. 15.321 a close-up of the dissipation functions with time at different Reynolds number 
is shown. A clear transition to turbulence is present, the appearance of progressively more 
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Figure 5.28: Isosurfaces of the squared current: side view. Two isosurfaces of the squared 
current at time t ~ 550 ta for a numerical simulation with 512x512x200 grid points and a 
Reynolds number TZ = 800 are shown. In partially transparent yellow is represented the 
isosurface at the value = 2.8 • 10^ while in red is the isosurface with = 8 • 10^, well 
below the value of the maximum of the squared current that at this time is j"^ = 3.6 • 10^. 
The red isosurface is always nested inside the yellow one, and in fact it appears as pink. 
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Figure 5.29: Isosurfaces of the squared current: top view. The same two isosurfaces of the 
previous figure are shown. In partiahy transparent yellow is represented the isosurface at 
the value = 2.8 • 10^ while in red is the isosurface with = 8 • 10^, weh below the value 
of the maximum of the squared current that at this time is = 3.6 • 10^. The isosurfaces 
are elongated along the axial direction, and the corresponding filling factor is small. 
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Figure 5.30: Left: Magnetic {Em) and kinetic {Ek) energies as a function of time for a long- 
time simulation with R = 400 and 256x256x100 grid points. After the initial reconnection 
event a statistically steady state is reached. Right: Ohmic (J) and viscous (O) dissipation 
as a function of time for a long-time simulation with R = 400 and 256x256x100 grid points. 
After the initial reconnection event a statistically steady state is reached. 

fine structure as the Reynolds number increases. Note that the dissipation is independent 
of the Reynolds number only beyond a threshold below which diffusion is important even 
at large scales, and consequently most of the energy is diffusively lost from the field at 
large scales before it can reach the small scales. In particular in the case of the sheared 
forcing considered in this simulation, if we decrease the Reynolds number below TZ = 50 
the diffusive time is so short that no instability can grow, and the system follows the linear 
saturation curves, showing no sign at all of nonlinear dynamics. 
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Figure 5.31: Ohmic dissipation (and the corresponding linear saturation curve) for three 
simulations performed with different Reynolds numbers TZ = 200, 400 and 800 . The first 
peak of the simulation with R = 800 reaches the value ~ 5100 (see Fig. I5.24p which goes 
beyond the range of the y axis. In the fully nonlinear stage ohmic dissipations roughly 
overlap for the three simulations, reaffirming our hypothesis that the dynamics is turbulent, 
so that beyond a threshold dissipation is independent of the Reynolds number. 



141 




Figure 5.32: Close-up of Ohmic dissipation as a function of time at different Reynolds 
numbers. A clear transition to turbulence is shown. In fact at higher Reynolds number 
correspond more fine structures becoming intermittent. 
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Chapter 6 

Nonlinear Analysis: Self-Organized 
Anisotropic Turbulence 

As already mentioned in § 14.31 one of the main original results of this work is 
the connection between the heating of closed magnetic structures in the solar corona and 
anisotropic turbulence. In an environment threaded by a strong axial magnetic field, Alfven 
wave packets traveling in opposite directions are naturally present. In coronal loops they 
are excited by photospheric motions on scales of the order of the cross-length of a granule 
(~ 1000 fcm). Their "collisions" give rise to a turbulent cascade which transfers energy 
from the large to the small scales, where it is finally dissipated, i.e. converted to heat and 
particle acceleration. Current sheets elongated in the axial direction are a feature found 
in most of the previous numerical simulations carried out to model a coronal loop. There 
have been many explanations about the formation of these structures. In the framework 
of anisotropic turbulence, these structures spontaneously develop. In fact the cascade of 
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energy mainly takes place in the orthogonal planes while it is strongly reduced in the axial 
direction, so the turbulent transfer of energy to smaller perpendicular scales results in 
sheet-like structures elongated in the direction of the main field. 



which we have given a review of the current state in § 12.2.41 As it usually happens in most 
of turbulence theory, the system studied is modeled as a three-dimensional three-periodic 
one. Periodicity leads to the cancellation of all the flux terms at the boundary surfaces. In 
a realistic model of a coronal loop, the axial direction cannot be modeled as periodic, and 
in particular the energy flux originating from the photospheric surfaces is just the source 
of energy for the system, which of course cannot be neglected. The inclusion of the energy 
flux terms at the boundaries gives to the system the property to be self- organized. In fact 
while the amplitude of the perturbations is usually imposed to be "small" on some physical 
ground, in a coronal loop it is the system itself that sets this amplitude, as described in the 
following sections. 

6.1 Self-Organization and Scalings 

In a turbulent system energy is injected into the system at the rate ej„ at the scale 
iin, it is subsequently scattered along the inertial range with the transfer rate e, and finally 
dissipated at the dissipation scale at the rate €d- For balance all these fluxes must be equal 



This equality still holds approximately when, as in our case, changes in time, because 
the more rapid dynamics at the small scales in the inertial and dissipation ranges adjust 



In this section we analyze our system using the theory of anisotropic turbulence, of 



Cm — e — 
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the spectrum rapidly compared to the slower dynamics of the large scales. In our case the 
value of ein is given by S (I4.93|) . the integral of the Poynting flux at the boundaries, and 
the transfer rate e is given by the anisotropic turbulence theory (see § I2.2.4|) . Both ej„ and 
e depend on the value of the internal fields, so that the balance equality (j6.ip fixes their 
value. 

In our numerical simulations the injection scale is ~ -^/4, that we approximate 
to i in the following calculations. Elsasser variables = u_l it b_L are the fundamental 
variables in turbulence theory. The nonlinear terms in the equations of reduced MHD ()2.53p - 
(j2.56p are symmetric in the exchange of these two variables 

(z^ • Vj_) z^ is symmetric for z"*" <^=^=- z~ (6-2) 

while the linear terms simply describe a wave propagation in opposite directions for the two 
fields z^. Furthermore, in our case, the boundary conditions (|2.6Up - (|2.6ip . that we rewrite 
for convenience, 

z~ + z+ = +2uq at z = (6.3) 

z+ + z~ = +2ul at z = L (6.4) 

are also symmetric in the exchange. In this way the rms values of z^ are approximately 
equal at all scales, i.e. indicating with 6z\ the rms value of the field z at the scale A 

^^A ~ ^^x ~ '^^A (6.5) 
This is equivalent to say that we expect the cross helicity 

= y d^x ux • = ^ j d^x (^\z+f - \z-f^ ~ (6.6) 
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to vanish. This is also confirmed by our linear analysis (see § (H) which has shown that the 
magnetic field is noticeably bigger than the velocity field, and by our numerical simulations 
(see Figure I5.1U|) . As the velocity field is smaller than the magnetic field we can also 
approximate in the Poynting flux integral (|4.93p the magnetic field at the injection scale d. 
with the Elsasser variable, i.e. 6h(^ ~ dzi, so that identifying the injections energy rate 
with S we have 



The energy injection rate is directly proportional to the rms of the Elsasser fields at the 
injection scale i. The velocity forcing at the boundaries injects energy with the rate (j6.7p 
at the scale ^i„ = £/4~£<<L smaller than the axial lenght L, developing a perturbation 

5zii « VA- 

The strongly reduced cascade along the axial direction implies that the wave pack- 
ets have size of the order L along this direction, while in the ortogonal plane a turbulent 
cascade with the characteristics described in § 12.2.41 takes place. We now apply the theory 
of anisotropic turbulence (see § I2.2.4p to have the value of the transfer energy rate e. 



We start from the anisotropic version of the Iroshnikov-Kraichnan theory described 



in equations (j2.134p - (j2.142p and characterized by the orthogonal spectrum E^^ oc fcj . 
Indicating with t\ = \/5z\ the eddy turn-over time and with r_4 ~ L/v^ the Alfvenic 
crossing time, the number of collisions for the fractional perturbation to build up to order 
unity is 



(6.7) 




(6.8) 
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so that the energy transfer time is given by 

Tx- NxTj^^tI/ta, (6.9) 

and the energy flux for unitary volume e is given by 

e~— -2 (6.10 

Usuafly the value of the spectral flux e or equivalently of the perturbation at the large 
scale are flxed to a "small" value on some physical ground. In our model this value is self- 
consistently established. Multiplying equation (I6.10p by the volume I'^L in order to have 
the total flux, and considering the injection scale \ = £ we have 

,^fL^J^—^—6zt (6.11) 

This is the flux of energy that leaves the injection scale towards the small scales, and that is 
supposed to have a constant value along all the inertial range. Equations {6. 7|) and h6.11\) 



show that the system is self-organized. In fact both and e depend by 6z£, the rms value 
of the Elsasser variables at the injection scale i, so that the internal dynamics depends by 
the injection of energy and the injection of energy depends by the internal dynamics. As 
shown in Figure \67\] the two fluxes satisfy the equality condition (|6.ip only for the value 



2 



(6.12) 

For a value of the perturbation smaller than (|6.12p the injection flux is higher than the 
dissipation one, so that the pertubation 5z£ grows toward the equilibrium value. On the 
opposite for higher values than (|6.12p the dissipation flux is higher than the injection one, 
decreasing the value of the perturbation 6ze. The system is dynamic and turbulent so that 
we will have fluctuations around this equilibrium value. 
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Figure 6.1: Injection energy rate Cm (I6.7p in dashed line and transfer energy rate along the 
inertial range e (|6.1U|) in continuous line, as a function of the rms value of the Elsasser vari- 
ables 6z£ at the injection scale £. 6zi can only fluctuate around the equilibrium value (j6.12p . 
The parameters used are £ = 1, L = 10, Up^ = 1 and vj( = 200, which have also been used 
to perform the numerical simulations. 



Substituting the equilibrium value (j6.12p in (|6.1ip or equivalenty (j6.7p . we obtain 
the average spectral flux 



As previously said (see ()6.ip ) the injection, transfer and dissipation rates balance, so that 
equation (|6.13p holds for all of them, in particular for the dissipation rate that is also 
the heating rate, i.e. the rate at which energy is converted to heat and particle acceleration 
at the dissipation scale £d- 

As shown in § 12.2.41 anisotropic turbulence is characterize by different regimes 
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to which correspond different slopes of the spectrum Ei.j_. The relative amplitude of the 
perturbation relatively to the axial Alfvenic velocity 6z£ jvj^ indicates in which regime is the 
system. In particular at smaller values of this ratio correspond higher slopes in the energy 
spectrum. In our case from equation (j6.12|) we have that 



\T I \ I ■ (6.14) 

VA \LJ \vaJ 

This means that changing the ratio v_A/uph, the relative amplitude of the perturbations 
changes, and we move among different regimes of anisotropic turbulence, in particular in- 
creasing the value of the axial magnetic field increaseas and the ratio (j6.14p decreases. 

Previous results can be generalized to the other scalings for anisotropic turbulence 
described in § 12.2.41 The different scalings are characterized by the different number of 
collisions A'^;^ that a wave packet must suffer for the perturbation to build up to order unity. 
Introducing the parameter a we can write 



Ta~A,t^~(^)""7A) (6.16) 



where t\ = X/5zx and = L/v_a so that for the energy transfer time we have 

L ) \5zx 
and finally for the transfer energy rate 

(A) (6.17) 

Tx \vaJ A" ^ ' 

Flux (I6.17P is supposed to be constant along the inertial range, so that we can consider its 
value at the injection scale \ = 
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(a) (b) 

Figure 6.2: (a) Injection energy rate Cin (|6.7p and transfer energy rates along the inertial 
range e (j6.18p for different values of the parameter a, as a function of the rms value of the 
Elsasser variables 5zi at the injection scale i. 5z£ can only fluctuate around the equilibrium 
value (j6.19p . (b) Coronal heating functions e* (j6.20p as a function of the axial Alfvenic 
velocity v_a for different values of the parameter a, showing that to higher values of a 
(corresponding to weaker turbulent regimes) coronal heating rates are more efficient. The 
parameters used both in (a) and (b) are £ = 1, L = 10, Uph = 1 and = 200. 

The injection energy rate is always given by (j6.7p . For a = 1 we have the anisotropic 
Kolmogorov spectrum , while for a = 2 we have the anisotropic IK spectrum kj_ and 
for a = 4 the SG94 spectrum based on 4-waves resonant interactions kj_ .In Figure [6. 2al 
energy fluxes e (j6.18p . for different values of the parameter a, and Cm (|6.7p are plotted as a 
function of the perturbation 6z£. There is always an equilibrium value where the two fluxes 
are equal, while for a higher value of the perturbation the transfer rate e is bigger than the 
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injection rate ej„ and for a lower value of the perturbation the transfer rate is smaher than 
the injection rate. From (I6.18P and (16. 7p we have that the equihbrium value for 6zi is given 
by 

a. 

(6-19) 

and substituting this value in (|6.18|) or (|6.7p we obtain the coronal heating scalings 

e*^f(-) vX^' (6.20) 



which are plotted, for different values of the parameter a, in Figure [6. 2b[ To higher values 
of the parameter a correspond higher heating rates. 

The regime of anisotropic turbulence is determined by the ratio 5zi/vjs,, which 
from ()6.19p is given by 

M.fn^f!V,A^. (6.21) 

As a > 1, the property that to higher values of the Alfvenic velocity corresponds a 
smaller relative perturbations is preserved, whatever the anisotropic turbulent regime is. 
This means that in general none of the different scalings (j6.20p corresponding to different 
values of a is valid for the whole range of the possible values of v^. In fact, starting form 
a low value of strong anisotropic turbulence will develop, corresponding to a spectral 
index k^^'^ and a coronal heating rate (j6.20p with a = 1 and e oc v^'^. Increasing vj^ the 
system transitions to an anisotropic IK regime characterized by a = 2, the spectral slope 

—9 5/3 / 

kj_ and e oc . Increasing still the value of the Alfvenic velocity will cause the system 
to transition in the regime characterized by a = 4, the spectral slope k^^"^ and e oc v^JJ' . 
So that the scalings for the coronal heating as a function of the Alfvenic velocity vj,, is 
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characterized by regions with increasing slopes in correspondence of increasing values of Vj[, 
each region characterized by the scaling ()6.20p with an increasing value of the parameter a. 



6.2 Timescales 

From equations (j6.15|) - (|6.2U|) we can derive the values of the eddy turn-over time 
Tx = X/6zx, the energy transfer rate Tx ()6.16p and the number of collisions at the scale A 
Nx (16.15p . After a few algebraic calculations we have 

, . _2_ _2a_ 

TX = (rxr^)^ (^-j , Tx = {t^ta)^ (^-j (6.22) 

and 

where t± = l/uph is the crossing time used as characteristic time to render the reduced 
MHD equations dimensionless, and Tj[ = L/v^x is the axial Alfvenic crossing time. For our 
system 

T± < TA (6.24) 

so that Nx ^ 1. The number of collisions has always the property to decrease at small 
scales, leading to a stronger turbulent regime. Except for a = 1 where correctly the energy 
transfer time and the eddy turn-over time are equal, for all the others regime (a > 2) 

Q — 1 Q Q — 1 

Tx I T±\ "+1 / A\ «+2 



> 1 (6.25) 



the energy transfer time is bigger than the eddy turn-over time, which is a characteristic of 
the Alfven effect at the basis of the anisotropic turbulence theory. 
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An interesting quantity is the energy transfer time at the large scales X = i which 
can be identified with the nonlinear timescale T^l, i.e. 

r^L = r^ = (r^T^)^ (6.26) 

1 

which in dimensionless units is T^l = ^2^^ • 

The nonlinear time Tjvl is also the time at which the linear stage transitions into 
the nonlinear stage. This results from the following calculations. The linear stage dynamics 
is described in § U and in particular we use the solution ()4.5ip - (|4.56p for the two-sided 
problem with velocity patterns and u^, respectively in the boundary planes z = and 
z = L. Averaging over time-scales bigger than the Alfvenic crossing time, the magnetic 
field for this solution can be approximated as 

(x, y, z,t)^— [u^ (x, y) - {x, y)] (6.27) 

Since during the linear stage the magnetic field has only large-scale components, we can 
approximate 

6ze ~ Uph — (6.28) 
and from (j4.87p we can write for the injection energy rate 

em ~ t^.A da (u^ - u°) • b_L ~ / da |u^ - u°|^ — ~ va uL — (6.29) 

J J TA TA 

From this equation we can also note the important property that during the linear stage the 
injection energy rate is always positive. The energy flux e along the inertial range is the rate 
at which the energy is leaving the large scales towards the small scales, and substituting 
the value (|6.28p for 6zi in (|6.18p we have for the energy transfer rate 
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Figure 6.3: Injection energy rate (16. 7p in continuous line and transfer energy rate along 
the inertial range e (j6.1U|) in dashed lines during the linear stage, as a function of time. The 
parameters used are £ = 1, L = W, Uph = 1 and vj\^ = 1000, to which corresponds a crossing 
Alfven time t_a = L/v_a = 10~^. 



a+2 



(6.30) 



Since a > 1 the injection of energy is higher than the dissipation, until time 



t = Tnl ~ (t" r^) -+i 



(6.31) 



where they become equal, and of course where the linear regime is no longer valid. 
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Chapter 7 

Conclusions and Discussion 

In this thesis I have presented the results of the numerical and analytical investi- 
gations of a reduced MHD realization of the Parker Scenario for coronal heating. This is 
the 3D extension of earlier 2D investigations by Einaudi, Velli, Politano & Pouquet [22j, 
Einaudi & VelU [21], Georgoulis, Velli & Einaudi t32j. 

Parker [68], [69], |70j . |71j was the first to propose that the X-ray corona is heated 
by dissipation at the many small current sheets forming in a coronal loop as a consequence 
of the continuous shuffling and intermixing of the footpoints of the field by the photospheric 
convection. Analogous to the linear analysis performed in §|l]the magnetic field results from 
a continuous mapping of the footpoints velocity pattern. He supposed that the magnetic 
field spontaneously produces tangential discontinuities (current sheets) which become in- 
creasingly severe with continued winding and interweaving, eventually producing intense 
magnetic dissipation by magnetic reconnection. This dissipation is characterized by bursts 
of rapid reconnection that he named "nanoflares" . 
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This picture is, however, "static" (as opposed to "dynamic" ) because the magnetic 
field-hnes are assumed to result from "quasi-static" displacements of the magnetic fields, 
and the current sheets are a consequence of the displacement of two neighboring field lines 
in opposite directions. The system does not show any dynamics, the field lines are pas- 
sively bent by the photospheric motions, and the formation of the current sheets is almost 
"geometric" . 

The first numerical simulation of the Parker Scenario was been performed in 1989 
by Mikic, Schnack &; Van Hoven [56], who found that current sheets formed, extended in 
the axial direction. Similar results were found later by Longcope & Sudan [53j and Hendrix 
&; Van Hoven [32] respectively in 1994 and 1996. The original "static" picture has remained 
the basis of all studies of the Parker Scenario. In particular, Longcope & Sudan [53] derive a 
scaling law for coronal heating as if this were due to magnetic reconnection of the large-scale 
field. Hendrix & Van Hoven ^2] also found current sheets extended in the axial direction 
and an inertial range but still thought that magnetic reconnection is responsible for the 
turbulence, describing it as "spontaneous unstably driven MHD turbulence" [32]. All these 
simulations used a relatively low resolution, i.e. a small Reynolds number, and could not 
observe the full turbulent dynamics. 

In this work, using high-resolution simulations, we have conclusively shown that 
the transfer of energy from the large-scales, where energy is continuously injected by pho- 
tospheric motions, towards the small scales is due to weak anisotropic turbulence. The 
small scale organizes into vortex-current sheets (a common feature in MHD turbulence) 
that result from the non-linear cascade. They may eventually break up due to tearing 
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and magnetic reconnection, but are certainly not primarily generated by a reconnection 
instability due to the large-scale field. Magnetic reconnection of these cascade-produced 
current sheets certainly plays a role in the dissipation of energy, i.e. in the conversion of 
the energy flowing along the inertial range into heat and particle acceleration. In any case 
the dissipation mechanism, which acts at the small scales, cannot be properly investigated 
with a fluid model because these are very well known to break at small scales and/or high 
frequencies. Thus, a kinetic investigation is appropriate. 

In my view the difficulty in realizing that the dynamics of a coronal loop is turbu- 
lent has been partially because a theory of weak turbulence for an MHD system embedded 
in a strong axial field has only recently been developed, mainly due to Sridhar & Goldreich 
(1994) [86] who studied this mechanism for the first time and find the spectrum. The 

debate followed with the remarks of Montgomery & Matthaeus (1995) |60j, while Goldreich 
& Sridhar (1997) [SI] found for the first time the spectrum and propose a multiple 
inertial range. Ng & Bhattacharjee (1997) [62], and Galtier et al. (2000,2002) [301 [M] then 
develop a kinetic model. 

Weak turbulence has been barely investigated numerically. Cho & Vishniac (2000), 
Maron & Goldreich (2001) and Miiller & Biskamp (2000) have performed three-dimensional 
simulations of strong MHD turbulence stretching the abilities of the fastest supercomputers. 
Simulations of weak turbulence are more difficult because each wavepacket must interact 
many times with oppositely directed ones before it cascades; by contrast, a wavepacket 
cascades in strong MHD turbulence in the time it takes to cross a single oppositely directed 
wavepacket. 



157 

Our investigations show that: 

The dynamics are described by weak turbulence. Because the system is embedded in 
a strong axial field the cascade takes place in the orthogonal planes and is strongly 
inhibited in the axial direction. Wave-packets are "long-lived" and they need to collide 
many times before transferring energy to smaller scales. 

The spectral slopes are in agreement with those derived from weak turbulence theory. 
The energy spectra slopes change from —2 for ~ 50 to ~ —3 for v^^ ~ 1000. The 
slope increase takes place while the rms of the orthogonal fields grows in absolute 
magnitude but its relative magnitude z'^/va decreases, in accordance with the theory. 

Small-scales structures do not homogeneously distribute in the planes but organize in 
vortex-current sheets extended in the axial direction where the energy flowing along 
the inertial range finally dissipates. 

In the framework of weak turbulence these extended structures arc naturally formed 
because the cascade in the axial direction is strongly inhibited. Hence no boundary 
layer is present. On the other hand in a real coronal loop a boundary layer, the 
transition region, is present. This result implies that it has others physical origins. 

Since the rate of injection of energy by photospheric motion depends not only on the 
photospheric velocity but also on the fields that develop into the computational box 
the system is self-organized. We have shown that the rms amplitude of the magnetic 
field at the large scales, which would grow linearly in time were it not for the non- 
linear dynamics, is fixed by the equality of the injection energy rate ei„ and the rate at 
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which energy flows from the large scales towards the small scales, which is the transfer 
energy rate along the inertial range e. 

• The weak turbulence and the self-organization of the system lead to the new coronal 
heating settlings ()6.20p 



where e is the rate of energy flowing from the injection scale towards the small scales 
along the inertial range where flnally dissipates. We have shown in ^ I6.1l that there is 
not a single scaling law, but that the dissipation rate depends on the regime of weak 
turbulence in which the system relaxes, described in the previous equation by the 
parameter a. We have also shown that to higher values of the Alfvenic velocity vji^ 
corresponds weaker turbulent regimes characterized by a higher value of the parameter 
a and thus leading to a more efficient heating rate. 

The system we have studied is a rather schematic model of a coronal loop. In 
particular we have not considered density or magnetic field variations along the axis of the 
loop. The least realistic assumptions are the boundaries, at which the energy is injected 
and hence are critical to model. Our model should be considered to apply to the extended 
central part of a loop, beyond the transition region, where the hypothesis of homogeneity 
along the axial direction is justified. 



ing velocity. A compressible code gives the possibility to use more realistic photospheric 
motions, including at the same time compressible effects into the dynamics. The first 
compressible simulations using an incompressible forcing velocity pattern (Dahlburg [13j) 





The first step towards a more complex model is to consider a more complex fore- 
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confirm the picture shown in this thesis. A future model for the numerical and theoretical 
investigation of a coronal loop dynamics should consider the underlying regions from the 
photosphere up to the transition region and of their mutual influences. 
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Part II 



Slow Solar Wind 
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Chapter 8 

The Slow Solar Wind Acceleration 

During the first part of my Pli.D. I fiave completed a work about tlie acceleration 
of the slow solar wind, which I had previously started for my "Laurea" thesis. This work 
has led to the publication of a paper (Rappazzo et al. [79]) which is reproduced in the 
Appendix. In the following a brief review is given. 

Although the association between the slow solar wind and the streamer belt (e.g. 
Gosling et al. [35j) and between the fast wind and the polar coronal holes is broadly rec- 
ognized, the mechanism which leads to the slow and fast acceleration is still a matter of 
debate. Einaudi et al. [21] developed an MHD model, with a current sheet embedded in 
a broader wake flow, that accounts for some of the typical features observed in the slow 
component of the solar wind. Reconnection of the magnetic field occurs at the current sheet 
and, in the non- linear regime, when the equilibrium magnetic field is substantially modi- 
fied, a Kelvin-Helmoltz instability develops, leading to the acceleration of density enhanced 
magnetic islands. 



164 

The solar streamer belt is a structure consisting of a magnetic configuration cen- 
tered on the current sheet, which extends above the cusp of a coronal helmet. The region 
underlying the cusp is made up of closed magnetic structures, with the cusp representing 
the point where separatrices between closed and open field lines intersect. Further from 
the Sun, at solar minimum, the streamer belt around the equator appears as a laminar 
configuration consisting of a thick plasma sheet with a density about 1 order of magnitude 
higher than the surrounding plasma, in which much narrower and complex structures are 
embedded. As a first approximation, moving from the center of the streamer in polar di- 
rections at radii greater than the radius of the cusp, the radial component of the magnetic 
field increases from zero, having opposite values on the two sides of the current sheet. As 
far as the flow distribution is concerned, the fast solar wind originates from the unipolar 
regions outside the streamer belt, while the slowest flows are located at the center of the 
sheet. 

One of the most interesting findings of the LASCO instrument onboard the Solar 
and Heliospheric Observatory (SOHO) spacecraft has been the observation of a continuous 
outflow of material in the solar streamer belt. An analysis performed using a difference 
image technique (Sheeley et al. [85], Wang et al. j9l]) has revealed the presence of plasma 
density enhancements, called blobs, accelerating away from the Sun. These plasmoids are 
seen to originate just beyond the cusps of helmet streamers as radially elongated structures 
a few percent denser than the surrounding plasma sheet, of approximately 1 Solar Radius 
(Rq) in length and 0.1 Rq in width. They are observed to accelerate radially outward 
maintaining constant angular spans at a nearly constant acceleration up to the velocity of 
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200-450 km/s, in the spatial region between about 5 and 30 Rq. It has been inferred that 
the blobs are tracers of the slow wind, being carried out by the ambient plasma flow. 

In my work [72] I have included in the previous model (Einaudi et al. [21] ) spherical 
geometry effects, taking into account either the radial divergence of the magnetic field 
lines and the average expansion suffered by a parcel of plasma propagating outward, using 
the Expanding Box Model (EBM), and the diamagnetic force due to the overall magnetic 
field radial gradients, the so-called melon-seed force. I have found that the values of the 
acceleration and density contrasts can be in good agreement with LASCO observations, 
provided the spherical divergence of the magnetic lines starts beyond a critical distance 
from the Sun and the initial stage of the formation and acceleration of the plasmoid is due 
to the cartesian evolution of MHD instabilities. This result provides a constraint on the 
topology of the magnetic field in the coronal streamer, which observationally is unknown. 
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Appendix 



In this appendix we have included the paper A. F. Rappazzo, M. Velh, G. Einaudi 
and R. B. Dahlburg "Diamagnetic and Expansion Effects on the Observable Properties of 
the Slow Solar Wind in a Coronal Streamer", ApJ, 2005, volume 633, part 1, pages 474-488 
http : //dx . doi . org/10 . 1086/431916. 
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